1 /***********************************************************************/
2 /* Open Visualization Data Explorer */
3 /* (C) Copyright IBM Corp. 1989,1999 */
4 /* ALL RIGHTS RESERVED */
5 /* This code licensed under the */
6 /* "IBM PUBLIC LICENSE - Open Visualization Data Explorer" */
7 /***********************************************************************/
8
9 #include <dxconfig.h>
10
11
12 /* SimplifySurface
13
14 is a module that approximates a triangulated surface and resamples data
15
16 */
17
18 /*********************************************************************/
19 /*
20
21 Author: Andre GUEZIEC, March/April 1997
22
23 */
24 /*********************************************************************/
25
26 #if defined(HAVE_STRING_H)
27 #include <string.h>
28 #endif
29
30 #if defined(HAVE_STRINGS_H)
31 #include <strings.h>
32 #endif
33
34 #include <dx/dx.h>
35 #include "simplesurf.h"
36
37 /*
38 this file contains the "guts" of the surface simplification procedure
39
40 As much as possible I am using simple types
41
42 For best memory management, I have converted my routines so that they
43 use the data explorer DXAllocate and DXFree calls
44
45 */
46
47 /*+--------------------------------------------------------------------------+
48 | |
49 | _dxfSimplifySurfaceDriver |
50 | |
51 +--------------------------------------------------------------------------+*/
52
53
_dxfSimplifySurfaceDriver(int nV,float * v,int nT,int * t,int dim_data,float * vertex_data,float tolerance,float data_tolerance,int * old_nV_after_manifold_conversion,int * old_nT_after_manifold_conversion,int * new_nV,float ** new_v,float ** new_vertex_data,int * new_nT,int ** new_t,int ** vertex_parents,int ** face_parents,int ** vertex_lut,int ** face_lut,float * old_positional_error,float ** new_positional_error,float ** face_normals,float ** old_face_areas,float ** new_face_areas,int preserve_volume,int simplify_boundary,int preserve_boundary_length)54 int _dxfSimplifySurfaceDriver(int nV, float *v,
55 int nT, int *t,
56 int dim_data,
57 float *vertex_data,
58 float tolerance,
59 float data_tolerance,
60 int *old_nV_after_manifold_conversion,
61 int *old_nT_after_manifold_conversion,
62 int *new_nV, float **new_v,
63 float **new_vertex_data,
64 int *new_nT, int **new_t,
65 int **vertex_parents,
66 int **face_parents,
67 int **vertex_lut, /* for fixing non-manifolds */
68 int **face_lut, /* same thing: for eliminating degenerate faces */
69 float *old_positional_error,
70 float **new_positional_error,
71 float **face_normals,
72 float **old_face_areas,
73 float **new_face_areas,
74 int preserve_volume,
75 int simplify_boundary,
76 int preserve_boundary_length)
77
78 /* all of the arrays that are ** must be allocated by the program
79 if any of such arrays is not allocated properly, than the
80 routine must return ERROR (zero). However, in case of an
81 error, there is no need to DXFree() such arrays, as they are
82 freed by the caller routine */
83
84 {
85 int success = 0;
86
87 int nE;
88
89 EdgeS *e = NULL;
90
91 Htable edge_table[1];
92
93 /* create a simplification data structure where all the relevant arrays
94 are copied */
95
96 SimpData simp_data[1];
97
98 bzero(simp_data, sizeof (SimpData));
99
100 bzero(edge_table, sizeof (Htable));
101
102 /* the first step is to convert the surface to a manifold, in case
103 Isosurface would have created non-manifold vertices (I know that happens sometimes)
104
105 It is also possible that a surface provided directly by the user would not
106 be a manifold */
107
108 if (_dxfToOrientedManifold(nV, v, nT, t, new_nV, new_v, new_nT, new_t, vertex_lut, face_lut,
109 &nE, edge_table, &e))
110 {
111 *old_nV_after_manifold_conversion = *new_nV;
112 *old_nT_after_manifold_conversion = *new_nT;
113
114 simp_data->nT = *new_nT; /* the number of vertices, edges and triangles */
115 simp_data->nE = nE; /* after conversion to manifold and before simplification */
116 simp_data->nV = *new_nV;
117 simp_data->data_dim = (vertex_data) ? dim_data: 0;
118
119 simp_data->preserve_volume = preserve_volume;
120 simp_data->simplify_boundary = simplify_boundary;
121 simp_data->preserve_boundary_length = preserve_boundary_length;
122
123 simp_data->tolerance = tolerance;
124 simp_data->data_tolerance = data_tolerance;
125
126 *old_face_areas = (float *) DXAllocateZero( simp_data->nT * sizeof (float));
127
128 *vertex_parents = (int *) DXAllocateZero( simp_data->nV * sizeof (int));
129
130 *face_parents = (int *) DXAllocateZero( simp_data->nT * sizeof (int));
131
132 /* part of the simplification data structure that references existing arrays */
133
134 simp_data->edge_array = e;
135 simp_data->e_table = edge_table;
136 simp_data->vert = (Vertex *) *new_v;
137 simp_data->triangles = *new_t;
138 simp_data->vertex_father = *vertex_parents;
139 simp_data->triangle_father = *face_parents;
140 simp_data->old_face_areas = *old_face_areas;
141 simp_data->vertex_lut = *vertex_lut;
142
143 if (*old_face_areas && *vertex_parents && *face_parents &&
144
145 (_dxfCreateSimplificationDataStructure(simp_data, vertex_data, old_positional_error)) &&
146
147 (_dxfSimplifyManifoldSurface(simp_data)))
148
149 {
150 /* free only the simplification data structures that are not exported by the driver routine,
151 to avoid running out of memory in _dxfBuildSimplifiedSurface(): */
152
153 _dxfPartialFreeSimplificationDataStructure(simp_data);
154
155 /* build a simplified surface using the vertex_father and triangle_father arrays */
156
157 if (_dxfBuildSimplifiedSurface(simp_data, new_nV, new_v, new_nT, new_t,
158 new_face_areas, face_normals,
159 new_positional_error, new_vertex_data))
160
161 success = 1;
162
163 _dxfFreeSimplificationDataStructure(simp_data, vertex_data, old_positional_error);
164
165 }
166
167
168 _dxfFreeHashTable(edge_table);
169 if (e) DXFree((Pointer)e);
170
171
172 }
173 /* if the conversion to oriented manifold fails, vertex_lut, new_v, new_t ... etc.
174 are freed in the caller routine */
175
176 return success;
177 }
178
179 /*+--------------------------------------------------------------------------+
180 | |
181 | _dxfToOrientedManifold |
182 | |
183 +--------------------------------------------------------------------------+*/
184
_dxfToOrientedManifold(int nV,float * v,int nT,int * t,int * new_nV,float ** new_v,int * new_nT,int ** new_t,int ** vlut,int ** flut,int * n_edges,Htable * e_table,EdgeS ** edges)185 int _dxfToOrientedManifold(int nV, float *v, int nT, int *t,
186 int *new_nV, float **new_v, int *new_nT, int **new_t, int **vlut, int **flut,
187 int *n_edges, Htable *e_table, EdgeS **edges)
188
189 {
190 int success = 0;
191
192 /* we first eliminate the degenerate triangles */
193
194 if (_dxfEliminateDegenerateTriangles(nV, nT, t, new_nT, new_t, flut))
195 {
196 /* after that the number of triangles does not change any more */
197
198 nT = *new_nT;
199
200 /* we then eliminate the standalone vertices */
201
202 if (_dxfEliminateStandaloneVertices(nV,v,new_nV,new_v,nT,*new_t,vlut))
203 {
204 /* we create a fathers array that is large enough to contain all
205 the triangle corners */
206
207 int *fathers = (int *) DXAllocateZero (3 * nT * sizeof (int));
208
209 if (fathers)
210 {
211 if (_dxfJoinTriangleCorners(nT, *new_t, e_table, n_edges, edges, fathers))
212 {
213
214 if (_dxfBuildOrientedManifold(new_nV, new_v, nT, *new_t, vlut, n_edges,
215 edges, e_table, fathers))
216
217 { success = 1;}
218 }
219
220 DXFree((Pointer)fathers);
221 }
222 }
223
224 }
225
226 return success;
227
228 }
229
230 /*+--------------------------------------------------------------------------+
231 | |
232 | _dxfJoinTriangleCorners |
233 | |
234 +--------------------------------------------------------------------------+*/
235
_dxfJoinTriangleCorners(int nT,int * t,Htable * e_table,int * n_edges,EdgeS ** edges,int * fathers)236 int _dxfJoinTriangleCorners(int nT, int *t, Htable *e_table, int *n_edges, EdgeS **edges, int *fathers)
237
238 /* we join two triangle corners only if they share a manifold edge
239
240 we define a manifold edge as an edge that has at most two incident triangles
241 and such that the two triangles are consistently oriented. Also, the two triangles
242 should not have the same three vertices */
243
244 {
245 int success = 0;
246
247 /* we first determine the edges that are manifold
248
249 in a second pass, we take all manifold edges and we join the corresponding triangle corners */
250
251
252 register int i;
253
254 register u_char j,k;
255
256 int nT3 = 3 * nT;
257
258 /* allocate an array big enough to contain all possible edges */
259
260 EdgeS *edge_array = (EdgeS *) DXAllocateZero (nT3 * sizeof (EdgeS));
261
262 if (edge_array)
263 {
264
265 EdgeS
266 the_edge[1],
267 *cur_edge,
268 *edge_in_table;
269
270 int
271 *incident_tri = SurfEdgeGetTriangles(the_edge),
272 *endpoints = SurfEdgeGetVertices (the_edge),
273 nE = 0,
274 total_nE = nT3,
275
276 /* initialize a pointer to the current triangle */
277
278 *the_triangle = t,
279 first_tri,
280 second_tri,
281 corner_v0first,
282 corner_v1first,
283 corner_v0second,
284 corner_v1second,
285 edge_paired_as_manifold = 0,
286 v0, /* surface edge end-points */
287 v1,
288 *triangles_in_table;
289
290 /* initialize the fathers array */
291
292 _dxfInitFather(nT3, fathers);
293
294 bzero((char *)e_table, sizeof(Htable)); /* reset the hash table */
295
296 /* the size of the table should be chosen such that 2**tbl_size ~ 2 * number of edges
297 for best performances. Here taking nT is plenty enough*/
298
299 if (!_dxfInitHashTable(e_table,
300 _dxfFlog2(nT),
301 (int (*)(KEYTYPE))_dxfEdgeSKey,
302 (int (*)(RECTYPE, RECTYPE))_dxfEdgeSFilter,
303 NULL))
304 goto error;
305
306
307 /* initialize the number of faces per edge to 1 */
308
309 SurfEdgeLabel(the_edge) = 1;
310
311 /* for each triangle, add the three edges to the hash table if they are not already
312 indexed in the hash table */
313
314 for (the_triangle = t, i=0; i<nT; i++, the_triangle += 3)
315 {
316
317 for (k=0;k<3;k++)
318
319 {
320 edge_paired_as_manifold = 0;
321 v0 = the_triangle[k]; /* surface edge end-points */
322 v1 = the_triangle[(k+1)%3];
323
324 if (v0>=v1)
325 {
326 endpoints[0] = v0; endpoints[1] = v1;
327 incident_tri[0] = i; incident_tri[1] = -1;
328 }
329 else
330 {
331 endpoints[0] = v1; endpoints[1] = v0;
332 incident_tri[0] = -1; incident_tri[1] = i;
333 }
334
335 edge_in_table = (EdgeS *)_dxfFindRecord(e_table,(KEYTYPE)the_edge);
336
337 /* if the edge is not already in the table we add it */
338
339 if (!edge_in_table)
340 {
341 /* copy this edge to the edge array */
342
343 memcpy(edge_array+nE, the_edge, sizeof(EdgeS));
344
345 /* add the edge to the hash table */
346
347 if (!_dxfAddRecord(e_table, (RECTYPE)(edge_array+nE), (KEYTYPE)(edge_array+nE)))
348 goto error;
349
350 /* increment the number of edges */
351
352 nE++;
353
354 }
355
356 else
357 { /* there is already an edge (v0,v1) in the table */
358
359 /* either one of the following two cases apply:
360
361 CASE ONE: we have already determined that the edge is not a manifold edge,
362
363 by giving it a label of zero. In that case,
364 we do not even add the edge to the list of edges,
365 but we put a label of zero in the edge slot.
366
367 and we increment the edge counter by one
368
369 */
370
371 if (SurfEdgeLabel(edge_in_table) == 0)
372 {
373
374 bzero(edge_array+nE, sizeof(EdgeS));
375
376 nE++;
377
378 }
379
380 /* CASE TWO: the edge is supposed to be a manifold edge with a label == 2
381
382 *but now* there is a third triangle in the picture.
383
384 so the edge cannot be manifold.
385 so we assign it a label of zero.
386
387 the problem is, we thought that we had a single edge, and we have
388 two, plus the new one is three. So we increment the edge counter by two.
389 */
390
391 else if (SurfEdgeLabel(edge_in_table) == 2)
392 {
393 /* again none of the three edges is a manifold edge */
394
395 SurfEdgeLabel(edge_in_table) = 0;
396
397 bzero(edge_array+nE, sizeof(EdgeS)); nE++;
398
399 bzero(edge_array+nE, sizeof(EdgeS)); nE++;
400
401 total_nE++; /*reinstate the edge that we thought was paired */
402 }
403
404 /* CASE THREE: we did not establish whether the edge is a manifold edge :
405
406 label == 1
407
408 but we have two triangles and we can establish it */
409
410 else if (SurfEdgeLabel(edge_in_table) == 1)
411 {
412 first_tri = i;
413
414 triangles_in_table = SurfEdgeGetTriangles(edge_in_table);
415
416 /* if the orientations of the faces do match,
417 we update the edge record to recognize that it is a valid edge,
418 and we will *later* join the corresponding vertices using the father array,
419 if it is established in the edge that the edge (v0,v1) *is indeed* manifold */
420
421 if ((v0>=v1) && triangles_in_table[0] == -1)
422 {
423
424 /* do the triangle orientations match ? */
425
426 /* so v0 and v1 are visited in that order in the_triangle.
427 they have to be visited in the opposite order in triangle_in_table */
428
429 second_tri = triangles_in_table[1];
430
431 /* we know that v1 appears first in second_tri (actually I am not sure what
432 this means because it is only true modulo circular permutation.
433
434 We can locate v1 in t1: */
435
436 j = 0; while (j<3 && t[second_tri*3+j] != v1) j++;
437
438 /* for the orientation to match, v0 must appear next, right after v1
439
440 also, we do not want the two triangles to have exactly the same three
441 vertices */
442
443 if (j < 3 && t[second_tri*3+(j+1)%3] == v0 &&
444 t[second_tri*3+(j+2)%3] != t[first_tri*3+(k+2)%3])
445 {
446 edge_paired_as_manifold = 1;
447 triangles_in_table[0] = first_tri;
448 }
449
450 }
451
452 else if ((v0<v1) && triangles_in_table[1] == -1)
453 {
454 second_tri = triangles_in_table[0];
455
456 j = 0; while (j<3 && t[second_tri*3+j] != v1) j++;
457
458 if (j < 3 && t[second_tri*3+(j+1)%3] == v0 &&
459 t[second_tri*3+(j+2)%3] != t[first_tri*3+(k+2)%3])
460 {
461 edge_paired_as_manifold = 1;
462 triangles_in_table[1] = first_tri;
463 }
464
465 }
466
467 if (edge_paired_as_manifold)
468 {
469 /* so it worked, and we indeed have a manifold edge (for the moment)
470
471 we set the edge label to 2 (manifold edge) in the table */
472
473 SurfEdgeLabel(edge_in_table) = 2;
474
475 total_nE--; /* substract the edge that is paired */
476 }
477 else
478 {
479 /* so there are (at least) two triangles sharing that edge, but they
480 are not matchable. so we know that the edge is not manifold
481
482 we set the edge label accordingly
483 */
484
485 SurfEdgeLabel(edge_in_table) = 0;
486
487 /* introduce a dummy edge in the edge array */
488
489 bzero(edge_array+nE, sizeof(EdgeS));
490
491 /* increment the number of edges */
492
493 nE ++;
494 }
495 } /* end of CASE 3 */
496 }
497 }
498 }
499
500 /* Now the manifold edges have been identified.
501 For each manifold edge, I join the triangle corners associated with this edge */
502
503 /* so I loop on the manifold edges */
504
505 for (cur_edge = edge_array, i=0 ; i<nE; i++, cur_edge++)
506 {
507 /* if the edge was established to be a manifold edge */
508
509 if (SurfEdgeLabel(cur_edge) == 2)
510 {
511
512 /* retrieve v0, v1, t0, t1 */
513
514 incident_tri = SurfEdgeGetTriangles(cur_edge);
515 endpoints = SurfEdgeGetVertices(cur_edge);
516
517
518 v0 = endpoints[0];
519
520 v1 = endpoints[1];
521
522 first_tri = incident_tri[0];
523 second_tri = incident_tri[1];
524
525 /* what are the corners associated with the first triangle incident_tri[0] ? */
526
527 /* ANSWER we loop on k until v0 is found */
528
529 k = 0 ; while (k<3 && t[first_tri*3+k] != v0) k++;
530
531 if (k < 3)
532 {
533 corner_v0first = first_tri * 3 + k;
534 corner_v1first = first_tri * 3 + (k+1)%3;
535
536 /* corners associated with the second tri ? */
537
538 /* ANSWER we loop on j until v1 is found */
539
540
541 j = 0; while (j<3 && t[second_tri*3+j] != v1) j++;
542
543 if (j < 3)
544 {
545 /* note that we have already verified that first_tri and second_tri
546 have different vertices when we established that the edge cur_edge was
547 manifold, so there is no need to check for this again */
548
549
550 corner_v1second = second_tri * 3 + j;
551 corner_v0second = second_tri * 3 + (j+1)%3;
552
553 /* join the two vertex pairs */
554 _dxfJoin(corner_v1second, corner_v1first, fathers);
555 _dxfJoin(corner_v0second, corner_v0first, fathers);
556 }
557 }
558 }
559 }
560
561 /* the hash table is not needed once the corners have been joined
562
563 because we are going to allocate many more arrays,
564 it is better to free the space now */
565
566 _dxfFreeHashTable(e_table);
567 bzero((char *)e_table, sizeof (Htable));
568
569 /* return the array of edges.
570 Edges need to be recomputed later in _dxfIndexEdges
571 we reallocate the array of edges at that time
572 */
573
574 *edges = edge_array;
575
576 /* return the number of edges */
577
578 *n_edges = nE;
579
580 success = 1;
581 }
582
583 return success;
584
585 error:
586
587 if (edge_array) DXFree((Pointer)edge_array);
588 _dxfFreeHashTable(e_table);
589
590 return 0;
591
592 }
593
594 /*+--------------------------------------------------------------------------+
595 | |
596 | _dxfBuildOrientedManifold |
597 | |
598 +--------------------------------------------------------------------------+*/
599
_dxfBuildOrientedManifold(int * nV,float ** v,int nT,int * t,int ** vlut,int * nE,EdgeS ** edges,Htable * e_table,int * fathers)600 int _dxfBuildOrientedManifold(int *nV, float **v, int nT, int *t, int **vlut,
601 int *nE, EdgeS **edges, Htable *e_table, int *fathers)
602 {
603
604 int
605 success = 0;
606
607 int
608 nT3 = 3 * nT,
609 new_nV = 0;
610
611 float
612 *new_v = NULL;
613
614 int
615 *new_t = NULL,
616 *old_vlut = *vlut,
617 *new_vlut = NULL;
618
619 register int i;
620
621 /* create a conversion table from old vertices to new vertices (triangle corners)
622 this is exactly the new set of triangles */
623
624
625 new_t = (int *) DXAllocateZero (nT3 * sizeof (int));
626
627 if (!new_t) goto error;
628
629 for (new_nV = i = 0; i < nT3; i++)
630 {
631 int corner_father;
632
633 /* perform path compression of the fathers at the same time */
634
635 if ((corner_father = _dxfFather(i,fathers)) == i)
636 {
637 new_t[i] = new_nV;
638 new_nV ++;
639 }
640
641 /* this works because the father always has lower index than the child */
642
643 else new_t[i] = new_t[corner_father];
644
645 /* corner_father should always be <= i */
646 }
647
648 /* new_nV determines the size of the new vertex array and new vlut array */
649
650 new_vlut = (int *) DXAllocateZero (new_nV * sizeof (int));
651
652 if (!new_vlut) goto error;
653
654 /* compute the new vertex look-up table */
655
656 /* loop through new_t. For each different vertex, look-up the old
657 vertex, and the look-up table of the old vertex */
658
659 for (i = 0; i < nT3; i++) if (i == fathers[i])
660
661 /* for each corner that is a corner group representative,
662 we have already assigned a new consecutive number.
663 We now update the entry in the new vertex look-up table by
664 retrieving the old vertex lookup entry of the vertex corresponding
665 to the corner */
666
667 new_vlut[new_t[i]] = old_vlut[t[i]];
668
669 /* we discard the old look-up table and replace the old look-up table
670 with the new lookup table */
671
672 if (*vlut)
673 {DXFree((Pointer)(*vlut));}
674
675 *vlut = new_vlut;
676
677 /* be careful, here, in case there is an error, new_vlut will be freed */
678
679
680 /* reshuffle the new vertex array */
681
682 new_v = (float *) DXAllocateZero (new_nV * sizeof (Vertex));
683
684 if (!new_v) goto error;
685
686 /* copy the vertices only if the corner is a father: only for each actual new vertex */
687
688 for (i = 0; i< nT3; i++)
689
690 if (i == fathers[i])
691
692 memcpy(new_v + 3 * new_t[i], *v + 3 * t[i], sizeof (Vertex));
693
694
695 /* discard the old vertex array and replace with the new vertex array */
696
697 if (*v) {DXFree((Pointer)(*v));}
698
699
700 *v = new_v;
701
702 /* be careful, here, in case there is an error, new_v will be freed
703 except that there cannot be an error after this point */
704
705 *nV = new_nV;
706
707
708 /* we overwrite the old triangle with the new triangles
709 and free the new triangles */
710
711 memcpy(t, new_t, nT3 * sizeof (int));
712
713 DXFree((Pointer) new_t);
714
715 new_t = NULL;
716
717 /* refill the edge hash table with the new vertex indices */
718
719 /* to start, reinitialize the hash table */
720 bzero((char *)e_table, sizeof (Htable));
721
722 if (_dxfIndexEdges(nE, edges, e_table, nT, t))
723 /* if this fails, we do not go to error, because there is nothing to
724 free (the hash table is freed in the error: of _dxfIndexEdges and so are the edges)
725 *v, *vlut ... will be freed later since they have been properly allocated
726 */
727 success = 1;
728
729 return success;
730
731 error:
732 if (new_t) {DXFree((Pointer)new_t);}
733
734 if (new_vlut)
735 {
736 DXFree((Pointer)new_vlut);
737 /* new_vlut was computed and vlut freed and replaced with new_vlut,
738 so we must set *vlut = NULL */
739 *vlut = NULL;}
740
741 if (new_v)
742 {
743 DXFree((Pointer)new_v);
744 /* new_v was computed and *v freed and replaced with new_v,
745 so we must set *v= NULL */
746 *v = NULL;
747 }
748
749 return 0;
750 }
751
752 /*+--------------------------------------------------------------------------+
753 | |
754 | _dxfIndexEdges |
755 | |
756 +--------------------------------------------------------------------------+*/
757
_dxfIndexEdges(int * n_edges,EdgeS ** new_edges,Htable * e_table,int nT,int * t)758 int _dxfIndexEdges(int *n_edges, EdgeS **new_edges, Htable *e_table, int nT, int *t)
759 {
760 /*
761 It is very important to count the edges again!
762
763 particularly, we may have implicitely joined edges that we thought were not joined */
764
765 EdgeS
766 *edges = NULL,
767 the_edge[1],
768 *edge_in_table;
769
770 int
771 nE = *n_edges,
772 *the_triangle = t,
773 *incident_tri = SurfEdgeGetTriangles(the_edge),
774 *endpoints = SurfEdgeGetVertices (the_edge),
775 second_tri,
776 v0, /* surface edge end-points */
777 v1,
778 *triangles_in_table;
779
780
781 register int i;
782
783 register u_char k;
784
785 /* this routine works for both a manifold and a non manifold surface
786 but the number of edges must be known in advance */
787
788 bzero((char *)e_table, sizeof(Htable)); /* reset the hash table */
789
790
791 /* initialize the number of faces per edge to 1 */
792
793 SurfEdgeLabel(the_edge) = 1;
794
795 /* reallocate the array of edges */
796
797 /* there may be fewer edges than nE, especially if some edges were
798 implicitely merged. However, there cannot be more than nE edges
799 So we are safe. However, it is *critical* to *not reallocate* after
800 the final hash_table is build in this procedure. Otherwise,
801 the edges may be relocated, and the hash-table entries will point
802 to NIL. */
803
804 *new_edges = (EdgeS *) DXReAllocate(*new_edges, nE * sizeof (EdgeS));
805
806 if (!(*new_edges)) goto error;
807
808 edges = *new_edges;
809
810 /* initialize the hash table with roughly as much entries as
811 the number of edges */
812
813 if (!_dxfInitHashTable(e_table,
814 _dxfFlog2(nE),
815 (int (*)(KEYTYPE))_dxfEdgeSKey,
816 (int (*)(RECTYPE, RECTYPE))_dxfEdgeSFilter,
817 NULL)) goto error;
818
819 for (the_triangle = t, nE = i = 0; i < nT; i++, the_triangle += 3)
820 {
821
822 for (k=0;k<3;k++)
823
824 {
825 v0 = the_triangle[k]; /* surface edge end-points */
826 v1 = the_triangle[(k+1)%3];
827
828
829 if (v0>=v1)
830 {
831 endpoints[0] = v0; endpoints[1] = v1;
832 incident_tri[0] = i; incident_tri[1] = -1;
833 second_tri = 0;
834 }
835 else
836 {
837 endpoints[0] = v1; endpoints[1] = v0;
838 incident_tri[0] = -1; incident_tri[1] = i;
839 second_tri = 1;
840 }
841
842 edge_in_table = (EdgeS *)_dxfFindRecord(e_table,(KEYTYPE)the_edge);
843
844 if (edge_in_table)
845 {
846
847 triangles_in_table = SurfEdgeGetTriangles(edge_in_table);
848
849 SurfEdgeLabel(edge_in_table) += 1;
850
851 triangles_in_table[second_tri] = i;
852
853 }
854 else
855 {
856
857 memcpy(edges+nE, the_edge, sizeof(EdgeS));
858
859
860 /* add the edge to the hash table */
861
862 if (!_dxfAddRecord(e_table, (RECTYPE)(edges+nE), (KEYTYPE)(edges+nE)))
863 goto error;
864
865 /* increment the number of edges */
866 nE++;
867
868 }
869 }
870 }
871
872 /* now the final number of edges is known.
873 I return this number to the caller routine */
874
875 *n_edges = nE;
876
877
878 return 1;
879
880 error:
881
882 _dxfFreeHashTable(e_table);
883
884 if (*new_edges) DXFree((Pointer)(*new_edges)); *new_edges = NULL;
885
886 return 0;
887 }
888
889 /*+--------------------------------------------------------------------------+
890 | |
891 | _dxfEliminateDegenerateTriangles |
892 | |
893 +--------------------------------------------------------------------------+*/
894
_dxfEliminateDegenerateTriangles(int nV,int nT,int * t,int * new_nT,int ** new_t,int ** flut)895 int _dxfEliminateDegenerateTriangles(int nV, int nT, int *t, int *new_nT, int **new_t, int **flut)
896 {
897 int success = 0; /* variable converted to one and returned if the procedure executes
898 successfully */
899
900 int
901 trisize = 3 * sizeof (int),
902 *t1 = t,
903 *new_tris = NULL,
904 *face_lut = NULL;
905
906 register int i;
907
908 /* initially, set the new number of triangles to zero */
909
910 *new_nT = 0;
911
912 /* 4/20/97: with respect to the original implementation, we now have the difference
913 that there may be a pre-existing "flut" face look-up table, which was created
914 for instance when quads were converted to triangles */
915
916
917 if (*flut == NULL)
918 {
919 face_lut = *flut = (int *) DXAllocateZero(nT * sizeof (int));
920
921 if (!(*flut)) goto error;
922
923 /* initialize the face_lut to point to the identity */
924
925 for (i=0;i<nT;i++) face_lut[i] = i;
926 }
927 else face_lut = *flut;
928
929
930 if (!(*flut)) goto error;
931
932 /* allocate more than enough new triangles and space for the new "flut":
933 we can't add any triangle above nT but we may eliminate some */
934
935 new_tris = *new_t = (int *) DXAllocateZero(nT * trisize);
936
937 if (new_tris && (*new_t))
938 {
939
940 /* take each old triangle in turn */
941
942 for (i=0 ; i<nT; i++, t1+=3)
943
944 /* if all three triangle indices are different, add the triangle to the new triangle list */
945
946 /* of course, all the indices should be 0 <= < nV
947
948 (we never know what sort of strange input people could create) */
949
950 if (t1[0] >= 0 && t1[0] < nV && t1[1] >= 0 && t1[1] < nV && t1[2] >= 0 && t1[2] < nV)
951
952 if (t1[0] != t1[1] && t1[1] != t1[2] && t1[0] != t1[2])
953 {
954 if (*new_nT < i)
955 face_lut[*new_nT] = face_lut[i]; /* i is always larger (>=) than *new_nT
956 so we will never overwrite a face_lut
957 entry that we will later need */
958
959 /* copy this triangle to the new list */
960
961 memcpy(new_tris + 3 * (*new_nT), t1, trisize);
962
963 *new_nT +=1;
964 }
965
966 /* now reallocate the new triangle array using the new number of triangles
967
968 reallocate also the face look up table */
969
970
971 if (nT != *new_nT)
972 {
973
974 *new_t = (int *) DXReAllocate(*new_t, *new_nT * trisize);
975
976 if (!(*new_t)) goto error;
977
978 *flut = (int *) DXReAllocate(*flut, *new_nT * sizeof (int));
979
980 if (!(*flut)) goto error;
981
982 }
983
984 /* don't forget to return "success" if there was no problem */
985
986 success = 1;
987 }
988
989 else goto error;
990
991
992 return success;
993
994 error:
995
996 if (*new_t) {DXFree((Pointer) (*new_t)); *new_t = NULL;}
997 if (*flut) {DXFree((Pointer) (*flut)); *flut = NULL;}
998
999 return 0;
1000 }
1001
1002 /*+--------------------------------------------------------------------------+
1003 | |
1004 | _dxfEliminateStandaloneVertices |
1005 | |
1006 +--------------------------------------------------------------------------+*/
1007
_dxfEliminateStandaloneVertices(int nV,float * v,int * new_nV,float ** new_v_exported,int nT,int * t,int ** vlut)1008 int _dxfEliminateStandaloneVertices(int nV, float *v, int *new_nV, float **new_v_exported,
1009 int nT, int *t, int **vlut)
1010 {
1011 /* the *new_v_exported notation stands for "Manifold set of Vertices" */
1012
1013 register int i, i1;
1014 register u_char k;
1015
1016 float *new_v = NULL;
1017
1018 int
1019 *new_tris = t,
1020 *vertex_lut = NULL,
1021 *old_to_new_vertex_index = NULL;
1022
1023 /* allocate the new arrays using the old counts since we potentially
1024 eliminate standalone vertices, the arrays cannot grow larger
1025 during this procedure (the vertex array may grow larger later on */
1026
1027 int *valence = (int *) DXAllocateZero(nV * sizeof (int));
1028
1029 if (!valence) goto error;
1030
1031 old_to_new_vertex_index = (int *) DXAllocateZero(nV * sizeof (int));
1032
1033 if (!old_to_new_vertex_index) goto error;
1034
1035 new_v = *new_v_exported = (float *) DXAllocateZero(3 * nV * sizeof (float));
1036
1037 if (!(*new_v_exported)) goto error;
1038
1039 vertex_lut = *vlut = (int *) DXAllocateZero(nV *sizeof (int));
1040
1041 if (!(*vlut)) goto error;
1042
1043 /* using the set of triangles, compute the valence of each vertex */
1044
1045 _dxfComputeVertexValence(nT, t, valence);
1046
1047 *new_nV = 0; /* reset the number of new vertices */
1048
1049 /* build the look-up table of vertices with non-zero valence: the
1050 remaining vertices */
1051
1052
1053 for (i=i1=0; i<nV; i++, i1+=3)
1054
1055 if (valence[i]>0)
1056 {
1057 vertex_lut[*new_nV] = i;
1058 old_to_new_vertex_index[i] = *new_nV;
1059
1060 /* copy the vertices that are not standalone in the new array */
1061
1062 memcpy(new_v + (*new_nV)*3, v + i1, 3 * sizeof (float));
1063
1064 *new_nV += 1; /* increment the count of new vertices */
1065 }
1066
1067 /* change the vertex indices in the triangle array */
1068
1069 for (i=i1=0; i<nT; i++, i1+=3)
1070
1071 /* for each triangle, change the index to its three vertices */
1072
1073 for (k=0; k<3; k++) new_tris[i1+k] = old_to_new_vertex_index[t[i1+k]];
1074
1075
1076 /* if there are fewer vertices than before, we reallocate the vertex array and the vertex lut */
1077
1078 if (*new_nV != nV)
1079 {
1080 *new_v_exported = (float *) DXReAllocate (*new_v_exported, *new_nV * sizeof (Vertex));
1081
1082 if (!(*new_v_exported)) goto error;
1083
1084 *vlut = (int *) DXReAllocate(*vlut, *new_nV * sizeof (int));
1085
1086 if (!(*vlut)) goto error;
1087 }
1088
1089 DXFree((Pointer) old_to_new_vertex_index);
1090 DXFree((Pointer) valence);
1091
1092 return 1;
1093
1094 error:
1095
1096 if (*new_v_exported) DXFree((Pointer) (*new_v_exported)); *new_v_exported = NULL;
1097 if (*vlut) DXFree((Pointer) (*vlut)); *vlut = NULL;
1098 if (valence) DXFree((Pointer) valence);
1099
1100 if (old_to_new_vertex_index) DXFree((Pointer) old_to_new_vertex_index);
1101
1102 return 0;
1103 }
1104
1105 /*+--------------------------------------------------------------------------+
1106 | |
1107 | dxfInitHashTable |
1108 | |
1109 +--------------------------------------------------------------------------+*/
1110
_dxfInitHashTable(Htable * table,int numbits,int (* fcode)(KEYTYPE),int (* filter)(RECTYPE,RECTYPE),KEYTYPE (* compute_key)(RECTYPE))1111 int _dxfInitHashTable(Htable *table, int numbits, int (*fcode)(KEYTYPE),
1112 int (*filter)( RECTYPE, RECTYPE), KEYTYPE (*compute_key)(RECTYPE))
1113 {
1114 register u_int i;
1115
1116 /* relation between size and number of bits */
1117 table->size = 1 << numbits;
1118
1119 /* the codes will be used as hash entries modulo table->mask */
1120 table->mask = table->size-1;
1121
1122 /* alloc the table of entries */
1123
1124 table->bucket = (Nodetype **) DXAllocateZero (sizeof(Nodetype *) * table->size);
1125
1126 if (!table->bucket) goto error;
1127
1128 for(i=0;i<table->size;i++)
1129 table->bucket[i] = (Nodetype *) NULL;
1130
1131 table->fcode = fcode; /* function that computes the code from the key */
1132 table->filter = filter;/* filtering function:
1133 returns 1 if the filter accepts the record
1134 returns 0 if the record is rejected, in which case the next
1135 record must be found and processed */
1136
1137
1138 if (compute_key == NULL) /* function that computes the key from the record */
1139 table->compute_key = _dxfIdentifyRecord2Key;
1140 else
1141 table->compute_key = compute_key;
1142
1143 /* initialize the first block of entries */
1144
1145 table->nodes_per_block = 1 << 10 ;
1146
1147 table->block = (Block *) NULL;
1148
1149 /* initialize the free node */
1150
1151 table->freenode = (Nodetype *) NULL;
1152
1153 return 1;
1154
1155 error:
1156
1157 bzero(table, sizeof (Htable));
1158 return 0;
1159 }
1160
1161 /*+--------------------------------------------------------------------------+
1162 | |
1163 | _dxfNewBlock |
1164 | |
1165 +--------------------------------------------------------------------------+*/
1166
_dxfNewBlock(Htable * table)1167 int _dxfNewBlock(Htable *table)
1168 {
1169 Block *newblock = NULL;
1170
1171 newblock = (Block *) DXAllocateZero (sizeof (Block));
1172
1173 if (!newblock) goto error;
1174
1175 newblock->next = table->block;
1176 newblock->rec = NULL;
1177
1178 table->block = newblock;
1179
1180 return 1;
1181
1182 error:
1183
1184 return 0;
1185 }
1186
1187 /*+--------------------------------------------------------------------------+
1188 | |
1189 | *_dxfgetnode |
1190 | |
1191 +--------------------------------------------------------------------------+*/
1192
_dxfgetnode(Htable * table)1193 Nodetype *_dxfgetnode(Htable *table)
1194 {
1195 Nodetype *retval;
1196
1197 register u_int i;
1198
1199 if (table->freenode)
1200 {
1201 retval = table->freenode;
1202 table->freenode = table->freenode->next;
1203 }
1204 else
1205 {
1206 /* create a new block */
1207
1208 /* this can potentially fail if there is not enough memory */
1209
1210 if (!_dxfNewBlock(table)) goto error;
1211
1212 /* allocate a new array of nodes */
1213
1214 table->block->rec= (RECTYPE) DXAllocateZero (sizeof(Nodetype)
1215 * table->nodes_per_block);
1216
1217 if (!table->block->rec) goto error;
1218
1219 /* link the nodes newly created */
1220
1221 for(retval= (Nodetype *) table->block->rec,
1222 i = 0; i < table->nodes_per_block-1; i++, retval++)
1223 {
1224 retval->next = retval+1 ;
1225 }
1226
1227 /* the last node has to point to NULL */
1228
1229 retval->next = (Nodetype *) NULL;
1230
1231 retval = (Nodetype *) table->block->rec ;
1232 table->freenode = retval+1;
1233 }
1234
1235 return(retval);
1236
1237 error:
1238
1239 return (NULL);
1240 }
1241
1242 /*+--------------------------------------------------------------------------+
1243 | |
1244 | *_dxfnewnode |
1245 | |
1246 +--------------------------------------------------------------------------+*/
1247
_dxfnewnode(Htable * table,RECTYPE rec,int code)1248 Nodetype *_dxfnewnode(Htable *table, RECTYPE rec, int code)
1249 {
1250
1251 Nodetype *retval = _dxfgetnode(table);
1252
1253 /* if _dxfgetnode(table) returns a NULL pointer, we should not try to access it */
1254
1255 if (retval)
1256 {
1257 retval->code = code;
1258 retval->rec = rec;
1259 }
1260 return(retval);
1261 }
1262
1263 /*+--------------------------------------------------------------------------+
1264 | |
1265 | _dxfHashTableReset |
1266 | |
1267 +--------------------------------------------------------------------------+*/
1268
_dxfHashTableReset(Htable * table)1269 void _dxfHashTableReset(Htable *table)
1270 {
1271 /* resets the hash table by cleaning all buckets and
1272 deleting all records */
1273
1274 Block
1275 *block=table->block,
1276 *next;
1277
1278 while(block)
1279 {
1280 next=block->next;
1281 if (block->rec) DXFree((Pointer) block->rec);
1282 DXFree((Pointer) block);
1283 block=next;
1284 }
1285 table->block = NULL;
1286
1287 /* clear buckets */
1288
1289 if (table->size > 0 && table->bucket)
1290 bzero((char *)table->bucket, table->size * sizeof(Nodetype *)) ;
1291
1292 return;
1293 }
1294
1295 /*+--------------------------------------------------------------------------+
1296 | |
1297 | _dxfAddRecord |
1298 | |
1299 +--------------------------------------------------------------------------+*/
1300
_dxfAddRecord(Htable * table,RECTYPE rec,KEYTYPE key)1301 int _dxfAddRecord(Htable *table, RECTYPE rec, KEYTYPE key)
1302 {
1303 Nodetype *node = NULL;
1304
1305 register int i;
1306
1307 int code= table->fcode(key);
1308
1309 node=_dxfnewnode(table,rec,code);
1310
1311 if (!node) goto error;
1312
1313 /* add the node containing the record to the appropriate bucket */
1314
1315 i=node->code & table->mask;
1316
1317 node->next = table->bucket[i];
1318 table->bucket[i]=node;
1319
1320 return 1;
1321
1322 error:
1323
1324 return 0;
1325 }
1326
1327 /*+--------------------------------------------------------------------------+
1328 | |
1329 | *_dxfFindRecord |
1330 | |
1331 +--------------------------------------------------------------------------+*/
1332
_dxfFindRecord(Htable * table,KEYTYPE key)1333 char *_dxfFindRecord(Htable *table, KEYTYPE key)
1334 {
1335 register int i;
1336
1337 table->search_code = table->fcode (key);
1338
1339 i = table->search_code & table->mask;
1340
1341 table->search_key = key;
1342
1343 table->search_node= table->bucket[i];
1344
1345 return (_dxfNextRecord(table));
1346 }
1347
1348 /*+--------------------------------------------------------------------------+
1349 | |
1350 | *_dxfNextRecord |
1351 | |
1352 +--------------------------------------------------------------------------+*/
1353
_dxfNextRecord(Htable * table)1354 char *_dxfNextRecord(Htable *table)
1355 {
1356 u_char found=0;
1357
1358 Nodetype *node=table->search_node;
1359
1360 /* if a filter function was provided, the first node with the same code AND such
1361 that the filter function comparing the node's record and the search key "agrees"
1362 is selected */
1363
1364 if (table->filter)
1365 {
1366 while (node && !found)
1367 {
1368 if (node->code == table->search_code &&
1369 table->filter(node->rec,table->compute_key(table->search_key)))
1370 found=1;
1371 else node = node->next;
1372 }
1373 }
1374
1375 /* if no filtering function was provided, the next node with the same code is
1376 retrieved */
1377
1378
1379 else
1380 while (node && !found)
1381 {
1382 if (node->code == table->search_code)
1383 found=1;
1384 else node = node->next;
1385 }
1386
1387 if (found)
1388 {
1389 table->search_node = node->next;
1390 return ((char *)node->rec);
1391 }
1392
1393 else return ((char *) NULL);
1394
1395 }
1396
1397 /*+--------------------------------------------------------------------------+
1398 | |
1399 | _dxfIdentifyRecord2Key |
1400 | |
1401 +--------------------------------------------------------------------------+*/
1402
1403 /* when a filtering function IS provided, we need to compare a current record with
1404 the search key. In the general case, a "compute_key" function must be provided
1405 to convert the key to a record or something that can be directly compared with
1406 a record
1407 in many cases the search key is the same as a record, or can be used without conversion.
1408 This is when the "dummy" function _dxfIdentifyRecord2Key is used.
1409 */
1410
_dxfIdentifyRecord2Key(RECTYPE rec)1411 KEYTYPE _dxfIdentifyRecord2Key(RECTYPE rec)
1412 {
1413 return (KEYTYPE) rec;
1414 }
1415
1416 /*+--------------------------------------------------------------------------+
1417 | |
1418 | _dxfFlog2 |
1419 | |
1420 +--------------------------------------------------------------------------+*/
1421 /*
1422 Log base 2 (integer form)
1423
1424 The following relation holds exactly
1425
1426 2 ^ Log2(N) <= N < 2 * 2 ^ Log2(N)
1427
1428 if N==0 then the result is by definition -1, but the true result
1429 is of course -infty
1430 */
1431
_dxfFlog2(int n)1432 int _dxfFlog2(int n)
1433 {
1434 int d=0;
1435 if (n==0) return (-1);
1436 while(n>1)
1437 {d++;
1438 n/=2;
1439 }
1440 return(d);
1441 }
1442
1443 /*+--------------------------------------------------------------------------+
1444 | |
1445 | _dxfEdgeSKey |
1446 | |
1447 +--------------------------------------------------------------------------+*/
1448
1449 /* compute the Hash-Table code from the key in the case of a surface edge */
1450
_dxfEdgeSKey(EdgeS * e)1451 int _dxfEdgeSKey(EdgeS *e)
1452 {
1453 return(e->v[1]^(e->v[0] << 8));
1454 }
1455
1456 /*+--------------------------------------------------------------------------+
1457 | |
1458 | _dxfEdgeSFilter |
1459 | |
1460 +--------------------------------------------------------------------------+*/
1461
_dxfEdgeSFilter(EdgeS * e1,EdgeS * e2)1462 int _dxfEdgeSFilter(EdgeS *e1, EdgeS *e2)
1463 {
1464 if ((e1->v[0] == e2->v[0]) &&
1465 (e1->v[1] == e2->v[1])) return(1);
1466 else
1467 return(0);
1468 }
1469
1470 /*+--------------------------------------------------------------------------+
1471 | |
1472 | _dxfEdgeSFilter |
1473 | |
1474 +--------------------------------------------------------------------------+*/
1475
_dxfFindEdge(Htable * e,int vA,int vB)1476 EdgeS *_dxfFindEdge(Htable *e, int vA, int vB)
1477 {
1478 EdgeS edge;
1479
1480 if (vA > vB)
1481 {
1482 edge.v[0]=vA;edge.v[1]=vB;
1483 }
1484 else
1485 {
1486 edge.v[0]=vB;edge.v[1]=vA;
1487 }
1488
1489 return (EdgeS *) _dxfFindRecord(e,(KEYTYPE) &edge);
1490
1491 }
1492
1493 /*+--------------------------------------------------------------------------+
1494 | |
1495 | _dxfInitArray |
1496 | |
1497 +--------------------------------------------------------------------------+*/
1498
_dxfInitArray(char * array,char * data,int n,int size)1499 void _dxfInitArray(char *array, char *data, int n, int size)
1500 {
1501
1502 register int i;
1503
1504 for(i=0;i<n;i++, array += size)
1505 memcpy(array, data, size);
1506
1507 return;
1508 }
1509
1510 /*+--------------------------------------------------------------------------+
1511 | |
1512 | _dxfInitFather |
1513 | |
1514 +--------------------------------------------------------------------------+*/
1515
_dxfInitFather(int n,int * father)1516 void _dxfInitFather(int n, int *father)
1517 {
1518 register int i;
1519
1520 for(i=0;i<n;i++) father[i] = i;
1521 }
1522
1523 /*+--------------------------------------------------------------------------+
1524 | |
1525 | _dxfFather |
1526 | |
1527 +--------------------------------------------------------------------------+*/
1528
_dxfFather(int i,int * father)1529 int _dxfFather(int i, int *father)
1530 {
1531 static int Ci, Fi;
1532
1533 /*--- path traversal ---*/
1534 for(Ci=i;(Fi=father[Ci])!=Ci;Ci=Fi);
1535
1536 /*--- path compression ---*/
1537 for(;(Fi=father[i])!=Ci;i=Fi) father[i] = Ci;
1538
1539 return(Ci);
1540 }
1541
1542 /*+--------------------------------------------------------------------------+
1543 | |
1544 | _dxfJoin |
1545 | |
1546 +--------------------------------------------------------------------------+*/
1547
_dxfJoin(int i,int j,int * father)1548 int _dxfJoin(int i, int j, int *father)
1549 {
1550 int rep_i, rep_j;
1551
1552 if((rep_i = _dxfFather(i,father)) != (rep_j = _dxfFather(j,father)))
1553 {
1554 /* the smallest representative becomes the father of
1555 the other representative */
1556
1557 if (rep_i < rep_j)
1558
1559 return (father[rep_j] = rep_i);
1560
1561 else if (rep_j < rep_i)
1562
1563 return (father[rep_i] = rep_j);
1564
1565 }
1566
1567 return -1;
1568 }
1569
1570 /*+--------------------------------------------------------------------------+
1571 | |
1572 | _dxfCreateSimplificationDataStructure |
1573 | |
1574 +--------------------------------------------------------------------------+*/
1575
_dxfCreateSimplificationDataStructure(SimpData * simp_data,float * vertex_data,float * old_positional_error)1576 int _dxfCreateSimplificationDataStructure(SimpData *simp_data, float *vertex_data, float *old_positional_error)
1577 {
1578 /* initialize all pointers except vx_data and err_volume to NULL */
1579
1580 simp_data->vx_data = vertex_data;
1581 simp_data->err_volume = old_positional_error;
1582
1583 simp_data->edge2index = NULL;
1584 simp_data->index2edge = NULL;
1585 simp_data->edge_father = NULL;
1586 simp_data->normal = NULL;
1587 simp_data->area = NULL;
1588 simp_data->compactness = NULL;
1589 simp_data->valence = NULL;
1590 simp_data->boundary_vert = NULL;
1591 simp_data->tol_volume = NULL;
1592 simp_data->vx_data_error = NULL;
1593 simp_data->vx_data_potential_values = NULL;
1594 simp_data->edge_heap = NULL;
1595
1596
1597 simp_data->edge2index = (int *) DXAllocateZero(simp_data->nE * sizeof (int));
1598
1599 if (!simp_data->edge2index) goto error;
1600
1601 simp_data->index2edge = (int *) DXAllocateZero(simp_data->nE * sizeof (int));
1602
1603 if (!simp_data->index2edge) goto error;
1604
1605 simp_data->edge_father = (int *) DXAllocateZero(simp_data->nE * sizeof (int));
1606
1607 if (!simp_data->edge_father) goto error;
1608
1609 /* initialize the father arrays */
1610
1611 _dxfInitFather(simp_data->nV, simp_data->vertex_father);
1612
1613 _dxfInitFather(simp_data->nE, simp_data->edge_father);
1614
1615 _dxfInitFather(simp_data->nT, simp_data->triangle_father);
1616
1617
1618 simp_data->normal = (Vertex *) DXAllocateZero(simp_data->nT * sizeof (Vertex));
1619
1620 if (!simp_data->normal) goto error;
1621
1622 simp_data->area = (float *) DXAllocateZero(simp_data->nT * sizeof (float));
1623
1624 if (!simp_data->area) goto error;
1625
1626 simp_data->compactness = (float *) DXAllocateZero(simp_data->nT * sizeof (float));
1627
1628 if (!simp_data->compactness) goto error;
1629
1630 /* compute the triangle normals, areas and compactness */
1631
1632 _dxfTrianglesNormalAreaCompactness(simp_data->nT, (Face *)simp_data->triangles, simp_data->vert,
1633 simp_data->normal, simp_data->area, simp_data->compactness);
1634
1635 /* save the areas in the "old_face_areas" array */
1636
1637 memcpy(simp_data->old_face_areas, simp_data->area, simp_data->nT * sizeof (float));
1638
1639
1640 /* use default values for the maximum angle between face normals
1641 and maximum ratio of compactness values */
1642
1643 simp_data->min_scalprod_nor = (float) cos ((double)SIMPLIFY_MAX_ANGLE_NOR);
1644
1645 simp_data->compactness_ratio = SIMPLIFY_COMPACTNESS_RATIO;
1646
1647
1648
1649 /* allocate and compute the vertex valence */
1650
1651 simp_data->valence = (u_short *) DXAllocateZero(simp_data -> nV * sizeof (u_short));
1652
1653 if (!simp_data->valence) goto error;
1654
1655 _dxfComputeVertexValence(simp_data->nT, (int *)simp_data->triangles, simp_data->valence);
1656
1657
1658 /* allocate and compute the "boundary_vert" array that flags boundary vertices with a flag value of 1 */
1659
1660 simp_data->boundary_vert = (u_char *) DXAllocateZero(simp_data -> nV * sizeof (u_char));
1661
1662 if (!simp_data->boundary_vert) goto error;
1663
1664 _dxfFlagBoundaryVertices(simp_data->nE, simp_data->edge_array, simp_data->boundary_vert);
1665
1666
1667 /* allocate and initialize the positional error */
1668
1669 /* there may or may not be positional error defined */
1670
1671 if (simp_data->err_volume)
1672 {
1673 /* there is positional error already defined */
1674
1675 float *old_error_volume = simp_data->err_volume;
1676
1677 simp_data->err_volume = (float *) DXAllocateZero(simp_data->nV * sizeof (float));
1678
1679 if (!simp_data->err_volume) goto error;
1680
1681 memcpy(simp_data->err_volume, old_error_volume, simp_data->nV * sizeof (float));
1682
1683 }
1684 else
1685 {
1686 simp_data->err_volume = (float *) DXAllocateZero(simp_data->nV * sizeof (float));
1687
1688 if (!simp_data->err_volume) goto error;
1689 }
1690
1691 simp_data->tol_volume = (float *) DXAllocateZero(simp_data->nV * sizeof (float));
1692
1693 if (!simp_data->tol_volume) goto error;
1694
1695 /* otherwise initialize the tolerance array */
1696
1697 _dxfInitArray((char *)simp_data->tol_volume,(char *) &simp_data->tolerance, simp_data->nV, sizeof (float));
1698
1699
1700 /* treat the vertex-related data */
1701
1702 /* there may or may not be vertex related data defined */
1703
1704 if (simp_data->vx_data && simp_data->data_dim > 0)
1705 {
1706 register int i;
1707
1708 float *old_vertex_data = simp_data->vx_data;
1709
1710 simp_data->vx_data = (float *) DXAllocateZero(simp_data->nV * simp_data->data_dim * sizeof (float));
1711
1712 if (!simp_data->vx_data) goto error;
1713
1714 simp_data->vx_data_error = (float *) DXAllocateZero(simp_data->nV * sizeof (float));
1715
1716 simp_data->vx_data_potential_values = (float *) DXAllocateZero(3 * simp_data->data_dim * sizeof (float));
1717
1718 if (!simp_data->vx_data_potential_values) goto error;
1719
1720 simp_data->vx_old_data = simp_data->vx_data_potential_values + simp_data->data_dim;
1721 simp_data->vx_new_data = simp_data->vx_old_data + simp_data->data_dim;
1722
1723 /* copy the vertex_data using the vertex_lut */
1724
1725 for (i=0;i<simp_data->nV;i++)
1726
1727 memcpy(simp_data->vx_data + i * simp_data->data_dim,
1728 old_vertex_data + simp_data->vertex_lut[i] * simp_data->data_dim,
1729 simp_data->data_dim* sizeof (float));
1730
1731 }
1732 else simp_data->vx_data = NULL;
1733
1734 /* allocate and initialize the edge heap */
1735
1736 simp_data->edge_heap = _dxfNewHeap(simp_data->nE);
1737
1738 if (!simp_data->edge_heap) goto error;
1739
1740 simp_data->get_lowest_weight_edge = _dxfRemoveLowestWeightEdgeFromHeap;
1741 simp_data->add_edge = _dxfAddEdge2Heap;
1742 simp_data->remove_edge = _dxfRemoveEdgeFromHeap;
1743
1744 simp_data->heap_initial_index = HeapInitialIndex(simp_data->edge_heap);
1745 simp_data->heap_outside_index = HeapOutsideIndex(simp_data->edge_heap);
1746
1747 /* initialize the heap indices */
1748
1749 _dxfInitArray((char *)simp_data->edge2index,(char *) &simp_data->heap_initial_index,
1750 simp_data->nE, sizeof (int));
1751
1752 /* initialize various counters: */
1753 simp_data->num_edg_collapsed = 0;
1754 simp_data->num_edg_weights = 0;
1755 simp_data->edg_rejected_4_topology = 0;
1756 simp_data->edg_rejected_4_geometry = 0;
1757 simp_data->edg_rejected_4_tolerance = 0;
1758
1759 /* initialize the last edge added to the heap: */
1760 simp_data->last_edge_added2heap = 0;
1761
1762 /* initialize the number of edges left to be tested: */
1763 simp_data->num_edg_remaining_4_testing = simp_data->nE;
1764
1765 simp_data->valence_max = SIMPLIFY_VALENCE_MAX;
1766
1767 return 1;
1768
1769 error:
1770
1771 if (simp_data->edge2index)
1772 {
1773 DXFree((Pointer) simp_data->edge2index);
1774 simp_data->edge2index = NULL;
1775 }
1776
1777 if (simp_data->index2edge)
1778 {
1779 DXFree((Pointer) simp_data->index2edge);
1780 simp_data->index2edge = NULL;
1781 }
1782
1783 if (simp_data->edge_father)
1784 {
1785 DXFree((Pointer) simp_data->edge_father);
1786 simp_data->edge_father = NULL;
1787 }
1788
1789 if (simp_data->normal)
1790 {
1791 DXFree((Pointer) simp_data->normal);
1792 simp_data->normal = NULL;
1793 }
1794
1795 if (simp_data->area)
1796 {
1797 DXFree((Pointer) simp_data->area);
1798 simp_data->area = NULL;
1799 }
1800
1801 if (simp_data->compactness)
1802 {
1803 DXFree((Pointer) simp_data->compactness);
1804 simp_data->compactness = NULL;
1805 }
1806
1807 if (simp_data->valence)
1808 {
1809 DXFree((Pointer) simp_data->valence);
1810 simp_data->valence = NULL;
1811 }
1812
1813 if (simp_data->boundary_vert)
1814 {
1815 DXFree((Pointer) simp_data->boundary_vert);
1816 simp_data->boundary_vert = NULL;
1817 }
1818
1819 if (simp_data->err_volume && simp_data->err_volume != old_positional_error)
1820 {
1821 DXFree((Pointer) simp_data->err_volume);
1822 simp_data->err_volume = old_positional_error;
1823 }
1824
1825 if (simp_data->tol_volume)
1826 {
1827 DXFree((Pointer) simp_data->tol_volume);
1828 simp_data->tol_volume = NULL;
1829 }
1830
1831 if (simp_data->vx_data && simp_data->vx_data != vertex_data)
1832 {
1833 DXFree((Pointer) simp_data->vx_data);
1834 simp_data->vx_data = vertex_data;
1835 }
1836
1837 if (simp_data->vx_data_error)
1838 {
1839 DXFree((Pointer) simp_data->vx_data_error);
1840 simp_data->vx_data_error = NULL;
1841 }
1842
1843 if (simp_data->vx_data_potential_values)
1844 {
1845 DXFree((Pointer)simp_data->vx_data_potential_values);
1846 simp_data->vx_data_potential_values = NULL;
1847 }
1848
1849 if (simp_data->edge_heap)
1850 {
1851 DXFree((Pointer) simp_data->edge_heap);
1852 simp_data->edge_heap = NULL;
1853 }
1854
1855 return 0;
1856 }
1857
1858 /*+--------------------------------------------------------------------------+
1859 | |
1860 | _dxfFreeSimplificationDataStructure |
1861 | |
1862 +--------------------------------------------------------------------------+*/
1863
_dxfFreeSimplificationDataStructure(SimpData * simp_data,float * vertex_data,float * old_positional_error)1864 int _dxfFreeSimplificationDataStructure(SimpData *simp_data, float *vertex_data, float *old_positional_error)
1865 {
1866
1867 if (simp_data->edge2index)
1868 {
1869 DXFree((Pointer) simp_data->edge2index);
1870 simp_data->edge2index = NULL;
1871 }
1872
1873 if (simp_data->index2edge)
1874 {
1875 DXFree((Pointer) simp_data->index2edge);
1876 simp_data->index2edge = NULL;
1877 }
1878
1879 if (simp_data->edge_father)
1880 {
1881 DXFree((Pointer) simp_data->edge_father);
1882 simp_data->edge_father = NULL;
1883 }
1884
1885 if (simp_data->normal)
1886 {
1887 DXFree((Pointer) simp_data->normal);
1888 simp_data->normal = NULL;
1889 }
1890
1891 if (simp_data->area)
1892 {
1893 DXFree((Pointer) simp_data->area);
1894 simp_data->area = NULL;
1895 }
1896
1897 if (simp_data->compactness)
1898 {
1899 DXFree((Pointer) simp_data->compactness);
1900 simp_data->compactness = NULL;
1901 }
1902
1903 if (simp_data->valence)
1904 {
1905 DXFree((Pointer) simp_data->valence);
1906 simp_data->valence = NULL;
1907 }
1908
1909 if (simp_data->boundary_vert)
1910 {
1911 DXFree((Pointer) simp_data->boundary_vert);
1912 simp_data->boundary_vert = NULL;
1913 }
1914
1915 if (simp_data->err_volume && simp_data->err_volume != old_positional_error)
1916 {
1917 DXFree((Pointer) simp_data->err_volume);
1918 simp_data->err_volume = old_positional_error;
1919 }
1920
1921 if (simp_data->tol_volume)
1922 {
1923 DXFree((Pointer) simp_data->tol_volume);
1924 simp_data->tol_volume = NULL;
1925 }
1926
1927 if (simp_data->vx_data && simp_data->vx_data != vertex_data)
1928 {
1929 DXFree((Pointer) simp_data->vx_data);
1930 simp_data->vx_data = vertex_data;
1931 }
1932
1933 if (simp_data->vx_data_error)
1934 {
1935 DXFree((Pointer) simp_data->vx_data_error);
1936 simp_data->vx_data_error = NULL;
1937 }
1938
1939 if (simp_data->vx_data_potential_values)
1940 {
1941 DXFree((Pointer)simp_data->vx_data_potential_values);
1942 simp_data->vx_data_potential_values = NULL;
1943 }
1944
1945 if (simp_data->edge_heap)
1946 {
1947 DXFree((Pointer) simp_data->edge_heap);
1948 simp_data->edge_heap = NULL;
1949 }
1950
1951 return 1;
1952 }
1953
1954 /*+--------------------------------------------------------------------------+
1955 | |
1956 | _dxfPartialFreeSimplificationDataStructure |
1957 | |
1958 +--------------------------------------------------------------------------+*/
1959
_dxfPartialFreeSimplificationDataStructure(SimpData * simp_data)1960 int _dxfPartialFreeSimplificationDataStructure(SimpData *simp_data)
1961 {
1962
1963 if (simp_data->edge2index)
1964 {
1965 DXFree((Pointer) simp_data->edge2index);
1966 simp_data->edge2index = NULL;
1967 }
1968
1969 if (simp_data->index2edge)
1970 {
1971 DXFree((Pointer) simp_data->index2edge);
1972 simp_data->index2edge = NULL;
1973 }
1974
1975 if (simp_data->edge_father)
1976 {
1977 DXFree((Pointer) simp_data->edge_father);
1978 simp_data->edge_father = NULL;
1979 }
1980
1981 if (simp_data->compactness)
1982 {
1983 DXFree((Pointer) simp_data->compactness);
1984 simp_data->compactness = NULL;
1985 }
1986
1987 if (simp_data->valence)
1988 {
1989 DXFree((Pointer) simp_data->valence);
1990 simp_data->valence = NULL;
1991 }
1992
1993 if (simp_data->boundary_vert)
1994 {
1995 DXFree((Pointer) simp_data->boundary_vert);
1996 simp_data->boundary_vert = NULL;
1997 }
1998
1999 if (simp_data->tol_volume)
2000 {
2001 DXFree((Pointer) simp_data->tol_volume);
2002 simp_data->tol_volume = NULL;
2003 }
2004
2005 if (simp_data->vx_data_error)
2006 {
2007 DXFree((Pointer) simp_data->vx_data_error);
2008 simp_data->vx_data_error = NULL;
2009 }
2010
2011 if (simp_data->vx_data_potential_values)
2012 {
2013 DXFree((Pointer)simp_data->vx_data_potential_values);
2014 simp_data->vx_data_potential_values = NULL;
2015 }
2016
2017 if (simp_data->edge_heap)
2018 {
2019 DXFree((Pointer) simp_data->edge_heap);
2020 simp_data->edge_heap = NULL;
2021 }
2022
2023 return 1;
2024 }
2025
2026 /*+--------------------------------------------------------------------------+
2027 | |
2028 | _dxfTrianglesNormalAreaCompactness |
2029 | |
2030 +--------------------------------------------------------------------------+*/
2031
_dxfTrianglesNormalAreaCompactness(int nT,Face * t,Vertex * v,Vertex * t_normal,float * t_area,float * t_compactness)2032 int _dxfTrianglesNormalAreaCompactness( int nT, Face *t, Vertex *v, Vertex *t_normal,
2033 float *t_area, float *t_compactness)
2034 {
2035 int success = 1;
2036
2037 static VertexD triVert[3];
2038 static double n[4];
2039
2040 register int i;
2041 register u_char j,k;
2042
2043 /* compute the face normals */
2044
2045 for(i=0;i<nT;i++,t++)
2046 {
2047 double
2048 sum_lengths_sq = 0.;
2049
2050 for (j=0;j<3;j++)
2051 for (k=0;k<3;k++)
2052 triVert[j][k] = v[t[0][j]][k];
2053
2054 _dxfTriangleNormalQR2D(triVert, n);
2055
2056 for (k=0;k<3;k++)
2057 t_normal[i][k] = n[k];
2058
2059 t_area[i] = n[3];
2060
2061 /* compute the squared lenghts of the three edges: */
2062
2063 for(j=0;j<3;j++)
2064 {
2065 double len_edg_sq = 0, tmp;
2066
2067 for(k=0;k<3;k++)
2068 {
2069 tmp = triVert[(j+1)%3][k] - triVert[j][k];
2070 len_edg_sq += tmp * tmp;
2071 }
2072
2073 sum_lengths_sq += len_edg_sq;
2074
2075 }
2076
2077 t_compactness[i] = FOUR_SQRT_3 * t_area[i] / sum_lengths_sq;
2078
2079 }
2080
2081 return success;
2082 }
2083
2084 /*+--------------------------------------------------------------------------+
2085 | |
2086 | _DxfTriangleNormalQR2D |
2087 | |
2088 +--------------------------------------------------------------------------+*/
2089
2090 /* computation of the triangle normal in double precision
2091
2092 n is an array
2093
2094 double n[4]
2095
2096 in n[3] the triangle area is returned
2097
2098
2099 this method is much more numerically stable than the cross product of the vectors.
2100
2101 it uses Householder orthogonalization of a set of two basis vectors chosen for the
2102 triangle. Experimental studies show that it is better to take both the shortest
2103 and longest edges of the triangle to form the basis
2104
2105 */
2106
2107 /* Turn off set but never used warnings for ("origin" below) */
2108 #ifdef sgi
2109 # pragma set woff 1552
2110 #endif
2111
_dxfTriangleNormalQR2D(VertexD * tri,double * n)2112 int _dxfTriangleNormalQR2D(VertexD *tri, double *n)
2113 {
2114
2115 static VertexD x1,x2;
2116
2117 int origin = 0;
2118
2119 _dxfTriangleBasisVectors(tri, x1, x2, origin);
2120
2121 return _dxfVectorProductQRD(x1,x2,n);
2122
2123 }
2124
2125 /*+--------------------------------------------------------------------------+
2126 | |
2127 | _dxfTriangleNormalQR2 |
2128 | |
2129 +--------------------------------------------------------------------------+*/
2130
2131 /* computation of the triangle normal in single precision
2132
2133 the triangle area is *not* computed */
2134
_dxfTriangleNormalQR2(Vertex * tri,Vertex n)2135 int _dxfTriangleNormalQR2(Vertex *tri, Vertex n)
2136 {
2137
2138 static float eps = 1e-8;
2139
2140 static Vertex x1,x2;
2141
2142 int origin = 0;
2143
2144 bzero((char *)n, 3 * sizeof (float));
2145
2146 _dxfTriangleBasisVectors(tri, x1, x2, origin);
2147
2148 {
2149
2150 float
2151 sign1 = (x1[0] >= 0.) ? 1. : -1.,
2152 sign2,
2153 norm1,
2154 norm2;
2155
2156 float beta1=0, beta;
2157
2158 static Vertex v1,v2;
2159
2160 register u_char k;
2161
2162 norm1 = NORM3(x1);
2163
2164 memcpy(v1,x1,sizeof(Vertex));
2165
2166 v1[0] += sign1 * norm1;
2167
2168 /* we have to record the sign of H1 x1: -sign1 */
2169
2170 if (norm1 >= eps)
2171
2172 /* we apply H1 to the vector x2 and replace x2 by H1 x2
2173 if the column is already zero, it is not necessary to use h1 */
2174 {
2175
2176 beta1 = - 2. / SCALPROD3(v1,v1);
2177
2178 beta = beta1 * SCALPROD3(v1,x2);
2179
2180 /* x2 = x2 + beta v1 */
2181
2182 for(k=0;k<3;k++)
2183 x2[k] += beta * v1[k];
2184 }
2185
2186 /* we create the v2 householder vector;
2187 it is by the way only a 2-vector */
2188
2189 {double norm2_d = x2[1]*x2[1] + x2[2] * x2[2];
2190
2191 if (norm2_d > 0.0) norm2 = (float) sqrt(norm2_d);
2192
2193 else norm2 = 0.0;}
2194
2195 /* we record the sign of H2 x2: -sign2 */
2196
2197 sign2 = (x2[1] >= 0.) ? 1. : -1.;
2198
2199 v2[0] = 0;
2200 v2[1] = x2[1] + sign2 * norm2;
2201 v2[2] = x2[2];
2202
2203
2204 if (norm2 >=eps)
2205
2206 {
2207 float beta2 = -2. *v2[2] / (v2[1]*v2[1] + v2[2]*v2[2]);
2208
2209 n[1] = beta2 * v2[1];
2210 n[2] = 1. + beta2 * v2[2];
2211 }
2212
2213
2214 if (norm1 >= eps)
2215 {
2216
2217 /* we apply H1 to the H2 e3 vector, and that's the normal
2218 up to a -1 factor:*/
2219
2220 beta = beta1 * SCALPROD3(v1, n);
2221
2222 /* n = n + beta v1 */
2223
2224 for(k=0;k<3;k++)
2225 n[k] += beta * v1[k];
2226 }
2227
2228
2229 if (sign1 * sign2 <0.)
2230
2231 for(k=0;k<3;k++)
2232
2233 n[k] = - n[k];
2234
2235 }
2236
2237
2238 return 1;
2239 }
2240
2241 #ifdef sgi
2242 # pragma reset woff 1552
2243 #endif
2244
2245 /*+--------------------------------------------------------------------------+
2246 | |
2247 | _dxfVectorProductQRD |
2248 | |
2249 +--------------------------------------------------------------------------+*/
2250
2251 /* compute the vector product using QR decomposition (which, to me, is the
2252 same as Householder orthogonalization. I guess the correct name is QR, since
2253 a Householder transformation represents a particular symmetry, and
2254 several Householder transformation must be combined to produce the final
2255 rotation + symmetry.
2256
2257 In the case of a 3 * 2 system, either one or two Householder transformations are required.
2258
2259 in n[3] store the area of the triangle defined by the origin
2260 and the two vectors x1 and x2 */
2261
2262
_dxfVectorProductQRD(VertexD x1,VertexD x2,double * n)2263 int _dxfVectorProductQRD(VertexD x1, VertexD x2, double *n)
2264 {
2265
2266 bzero((char *)n, 4 * sizeof (double));
2267
2268 {
2269 static double eps = 1e-15;
2270
2271 /* do a QR decomposition of the matrix (x1,x2) */
2272
2273 /* the first householder vector x1 is equal to :
2274 x1 + sign(x11) |x1| (1 0 0)^t
2275 we will call H1 the corresponding Householder reflexion
2276 */
2277
2278
2279 double
2280 sign1 = (x1[0] >= 0.) ? 1. : -1.,
2281 sign2,
2282 norm1,
2283 norm2;
2284
2285 double beta1=0, beta;
2286
2287 static VertexD v1,v2;
2288
2289 register u_char k;
2290
2291 norm1 = NORM3(x1);
2292
2293 memcpy(v1,x1,sizeof(VertexD));
2294
2295 v1[0] += sign1 * norm1;
2296
2297 /* we have to record the sign of H1 x1: -sign1 */
2298
2299 if (norm1 >= eps)
2300
2301 /* we apply H1 to the vector x2 and replace x2 by H1 x2
2302 if the column is already zero, it is not necessary to use h1 */
2303 {
2304 /* beta = -2/ (v1^t v1)* (v1^x2) */
2305
2306
2307 beta1 = - 2. / SCALPROD3(v1,v1);
2308
2309 beta = beta1 * SCALPROD3(v1,x2);
2310
2311 /* x2 = x2 + beta v1 */
2312
2313 for(k=0;k<3;k++)
2314 x2[k] += beta * v1[k];
2315 }
2316
2317 /* we create the v2 householder vector;
2318 it is by the way only a 2-vector */
2319
2320 norm2 = x2[1]*x2[1] + x2[2] * x2[2];
2321
2322 if (norm2 > 0.0) norm2 = sqrt(norm2);
2323
2324 /* we record the sign of H2 x2: -sign2 */
2325
2326 sign2 = (x2[1] >= 0.) ? 1. : -1.;
2327
2328 v2[0] = 0;
2329 v2[1] = x2[1] + sign2 * norm2;
2330 v2[2] = x2[2];
2331
2332 /* compute the area of the triangle */
2333
2334 n[3] = norm1 * norm2 /2.;
2335
2336
2337 if (norm2 >=eps)
2338
2339 /* we apply H2 to the e3 vector and put the result directly in n */
2340 {
2341 /* beta2 = -2 /v2^t v2 * v2^t e3 */
2342 double beta2 = -2. *v2[2] / (v2[1]*v2[1] + v2[2]*v2[2]);
2343
2344 n[1] = beta2 * v2[1];
2345 n[2] = 1. + beta2 * v2[2];
2346 }
2347
2348
2349 if (norm1 >= eps)
2350 {
2351
2352 /* we apply H1 to the H2 e3 vector, and that's the normal
2353 up to a -1 factor:*/
2354
2355 beta = beta1 * SCALPROD3(v1, n);
2356
2357 /* n = n + beta v1 */
2358
2359 for(k=0;k<3;k++)
2360 n[k] += beta * v1[k];
2361 }
2362
2363
2364 /* There is a unique QR decomposition such that the diagonal of R
2365 is made of positive elements and Q is a rotation.
2366 The last vector of Q is going to be n.
2367 if ( R(1,1)*R(2,2) < 0 ) we have to invert n */
2368
2369 if (sign1 * sign2 <0.)
2370
2371 for(k=0;k<3;k++)
2372
2373 n[k] = - n[k];
2374
2375 }
2376
2377 return 1;
2378
2379
2380 }
2381
2382 /*+--------------------------------------------------------------------------+
2383 | |
2384 | _dxfFlagBoundaryVertices |
2385 | |
2386 +--------------------------------------------------------------------------+*/
2387
2388 /* as the name indicates, mark the boundary vertices with a
2389 value of 1 in the array "boundary_vert"
2390 the interior vertices are marked with 0 */
2391
_dxfFlagBoundaryVertices(int nE,EdgeS * edges,u_char * boundary_vert)2392 int _dxfFlagBoundaryVertices(int nE, EdgeS *edges, u_char *boundary_vert)
2393 {
2394 /* loop through edges.
2395
2396 if an edge is a boundary edge, then mark both vertices as boundary vertices */
2397
2398 register int i;
2399
2400 for (i=0; i<nE; i++, edges++)
2401
2402
2403 if (SurfEdgeLabel(edges) == 1)
2404 {
2405 int *endpoints = SurfEdgeGetVertices(edges);
2406
2407 boundary_vert[endpoints[0]] = boundary_vert[endpoints[1]] = 1;
2408
2409 }
2410
2411 return 1;
2412 }
2413
2414 /*
2415
2416 static functions for Heap processing
2417
2418 */
2419
2420 static void HeapSwitch(int i, int j, Heap *h);
2421 static void HeapDown(int father, Heap *h);
2422 static void HeapUp(Heap *h);
2423 static int HeapDel(int i, Heap *h);
2424
2425 /*+--------------------------------------------------------------------------+
2426 | |
2427 | HeapSwitch |
2428 | |
2429 +--------------------------------------------------------------------------+*/
2430
HeapSwitch(int i,int j,Heap * h)2431 static void HeapSwitch(int i, int j, Heap *h)
2432 {
2433 /* switch two elements of the Heap, and update the permutation
2434 and inverse permutation arrays */
2435
2436 if(i!=j && i<h->n && j<h->n)
2437 {
2438 float key;
2439 int pos;
2440
2441 key = h->key[i]; h->key[i] = h->key[j]; h->key[j] = key;
2442 pos = h->invp[i]; h->invp[i] = h->invp[j]; h->invp[j] = pos;
2443
2444 h->perm[h->invp[i]] = i;
2445 h->perm[h->invp[j]] = j;
2446 }
2447 }
2448
2449
2450 /*+--------------------------------------------------------------------------+
2451 | |
2452 | HeapDown |
2453 | |
2454 +--------------------------------------------------------------------------+*/
2455
HeapDown(int father,Heap * h)2456 static void HeapDown(int father, Heap *h)
2457 {
2458
2459 /* trickle an element down the Heap until
2460 the Heap condition is verified
2461
2462 used when deleting an element from the Heap
2463 */
2464
2465
2466 if(0<=father && father<h->last-1)
2467 {
2468 int left,right,child;
2469
2470 while((left=LeftSon(father))<h->last)
2471 {
2472 child = father;
2473 if(h->key[left]<h->key[child]) child = left;
2474 if((right=RightSon(father))<h->last &&
2475 h->key[right]<h->key[child]) child = right;
2476 if(child==father)
2477 break;
2478 HeapSwitch(father,child,h);
2479 father = child;
2480 }
2481 }
2482 }
2483
2484 /*+--------------------------------------------------------------------------+
2485 | |
2486 | HeapUp |
2487 | |
2488 +--------------------------------------------------------------------------+*/
2489
HeapUp(Heap * h)2490 static void HeapUp(Heap *h)
2491 {
2492
2493 /* trickle the last heap element up the Heap until the heap
2494 condition is verified
2495
2496 used when adding an element to the heap */
2497
2498 int child,father;
2499
2500 child = h->last-1;
2501
2502 while(child>0 && h->key[(father=(child-1)/2)]>=h->key[child])
2503 {
2504 HeapSwitch(father,child,h);
2505 child = father;
2506 }
2507 }
2508
2509 /*+--------------------------------------------------------------------------+
2510 | |
2511 | HeapDel |
2512 | |
2513 +--------------------------------------------------------------------------+*/
2514
HeapDel(int i,Heap * h)2515 static int HeapDel(int i, Heap *h)
2516 {
2517
2518 /* delete an element inside a Heap using its actual Heap index
2519 this procedure is not meant to be called outside this file */
2520
2521 if (0<=i && i<h->last)
2522 {
2523 h->last--;
2524 HeapSwitch(i,h->last,h);
2525 HeapDown(i,h);
2526 return 1;
2527 }
2528 else return 0;
2529 }
2530
2531 /*+--------------------------------------------------------------------------+
2532 | |
2533 | _dxfHeapDelete |
2534 | |
2535 +--------------------------------------------------------------------------+*/
2536
_dxfHeapDelete(int i,Heap * h)2537 int _dxfHeapDelete(int i, Heap *h)
2538 {
2539 /* delete a Heap element using the original index it was assigned
2540 when first added to the Heap.
2541 This is the preferred procedure for deleting Heap elements
2542 (together with HeapDelMin())*/
2543
2544 if(0<=i && i<h->n)
2545 return HeapDel(h->perm[i],h);
2546
2547 else return 0;
2548 }
2549
2550 /*+--------------------------------------------------------------------------+
2551 | |
2552 | _dxfHeapDelMin |
2553 | |
2554 +--------------------------------------------------------------------------+*/
2555
_dxfHeapDelMin(Heap * h)2556 int _dxfHeapDelMin(Heap *h)
2557 {
2558
2559 /* retrieve the Heap element with the smallest key
2560 (first element in the Heap), take it off the Heap,
2561 and reorganize the Heap */
2562
2563 if(h->last>0)
2564 {
2565 HeapDel(0,h);
2566 return(h->invp[h->last]);
2567 }
2568 else
2569 return(-1); /* heap is empty */
2570 }
2571
2572 /*+--------------------------------------------------------------------------+
2573 | |
2574 | *_dxfNewHeap |
2575 | |
2576 +--------------------------------------------------------------------------+*/
2577
_dxfNewHeap(int n)2578 Heap *_dxfNewHeap(int n)
2579 {
2580
2581 /* create a new Heap: allocate and initialize */
2582
2583 Heap *heap = (Heap*)0;
2584
2585 if(n>0)
2586 {
2587 unsigned long bytes = _dxfHeapSize(n);
2588 char *block = (char *) DXAllocateZero(bytes);
2589
2590 if (!block) goto error;
2591
2592 heap = (Heap*)block;
2593 heap->n = n;
2594 heap->key = (float*)(block+sizeof(Heap));
2595 heap->perm = (int*)(block+sizeof(Heap)+n*sizeof(float));
2596 heap->invp = (int*)(block+sizeof(Heap)+n*sizeof(float)+n*sizeof(int));
2597
2598 _dxfResetHeap(heap);
2599 }
2600
2601 return(heap);
2602
2603 error:
2604
2605 return NULL;
2606 }
2607
2608 /*+--------------------------------------------------------------------------+
2609 | |
2610 | _dxfHeapSize |
2611 | |
2612 +--------------------------------------------------------------------------+*/
2613
_dxfHeapSize(int n)2614 unsigned long _dxfHeapSize(int n)
2615 {
2616 /* used by NewHeap(), compute the size that is needed
2617 given the maximum number of elements that can be contained
2618 by the Heap */
2619
2620
2621 unsigned long bytes = sizeof(Heap)+n*sizeof(float)+2*n*sizeof(int);
2622
2623 return(bytes);
2624 }
2625
2626 /*+--------------------------------------------------------------------------+
2627 | |
2628 | _dxfResetHeap |
2629 | |
2630 +--------------------------------------------------------------------------+*/
2631
_dxfResetHeap(Heap * heap)2632 void _dxfResetHeap(Heap *heap)
2633 {
2634
2635 /* without freeing the memory, re-initialize a Heap for
2636 subsequent use */
2637
2638
2639 heap->last = 0;
2640
2641 if(heap->n>0)
2642 {
2643 register int i;
2644
2645 for(i=0;i<heap->n;i++)
2646 {
2647 heap->key[i]=0.0;
2648 heap->perm[i] = heap->invp[i] = i;
2649 }
2650 }
2651
2652 return;
2653 }
2654
2655 /*+--------------------------------------------------------------------------+
2656 | |
2657 | _dxfHeapAdd |
2658 | |
2659 +--------------------------------------------------------------------------+*/
2660
_dxfHeapAdd(float key,Heap * h)2661 int _dxfHeapAdd(float key, Heap *h)
2662 {
2663 int indx = -1;
2664
2665 if(h->last<h->n) /* don't want to enlarge the heap if h->last==h->n ? */
2666 {
2667 int last = h->last++;
2668
2669 indx = h->invp[last];
2670 h->key[last] = key;
2671 HeapUp(h);
2672 }
2673
2674 return(indx);
2675 }
2676
2677 /*+--------------------------------------------------------------------------+
2678 | |
2679 | _dxfRemoveEdgeFromHeap |
2680 | |
2681 +--------------------------------------------------------------------------+*/
2682
2683
2684 /* remove a surface edge from the Heap using the edge's original
2685 index */
2686
2687
_dxfRemoveEdgeFromHeap(SimpData * simp_data,int index,int the_edge)2688 int _dxfRemoveEdgeFromHeap(SimpData *simp_data, int index, int the_edge)
2689 {
2690 return _dxfHeapDelete(index, simp_data->edge_heap);
2691 }
2692
2693 /*+--------------------------------------------------------------------------+
2694 | |
2695 | _dxfAddEdge2Heap |
2696 | |
2697 +--------------------------------------------------------------------------+*/
2698
2699 /* add a surface edge to the heap and record the index associated
2700 with the edge */
2701
_dxfAddEdge2Heap(SimpData * simp_data,int the_edge)2702 int _dxfAddEdge2Heap(SimpData *simp_data, int the_edge)
2703 {
2704 int
2705 added = 0;
2706
2707
2708 /* verify that the Heap is not already full and that
2709 the edge does not have an index that forbids adding to the Heap */
2710
2711
2712 if (!HeapFull(simp_data->edge_heap) && simp_data->edge2index[the_edge] >=0)
2713 {
2714 /* compute the representatives of the edge vertices */
2715 int
2716 *v = SurfEdgeGetVertices(simp_data->edge_array+the_edge),
2717 v0 = FATHER(v[0], simp_data->vertex_father),
2718 v1 = FATHER(v[1], simp_data->vertex_father);
2719
2720 /* 1) to compute the key to the heap, we choose the edge length to start with */
2721 float
2722 weight = DIST3(simp_data->vert[v0], simp_data->vert[v1]);
2723
2724 /* 1.1) we add the sum of the errors at the two vertices */
2725
2726 weight += simp_data->err_volume[v0] + simp_data->err_volume[v1];
2727
2728 /* 2) add the edge to the heap */
2729
2730
2731 simp_data->index2edge[
2732 simp_data->edge2index[the_edge]=
2733 _dxfHeapAdd(
2734 weight,
2735 simp_data->edge_heap)
2736 ] = the_edge;
2737 added = 1;
2738
2739 }
2740
2741 return added;
2742
2743 }
2744
2745 /*+--------------------------------------------------------------------------+
2746 | |
2747 | _dxfRemoveLowestWeightEdgeFromHeap |
2748 | |
2749 +--------------------------------------------------------------------------+*/
2750
_dxfRemoveLowestWeightEdgeFromHeap(SimpData * simp_data)2751 int _dxfRemoveLowestWeightEdgeFromHeap(SimpData *simp_data)
2752 {
2753 int
2754 idx = _dxfHeapDelMin(simp_data->edge_heap);
2755
2756 if (idx != -1)
2757
2758 return simp_data->index2edge[idx];
2759
2760 else
2761 return -1;
2762 }
2763
2764 /*+--------------------------------------------------------------------------+
2765 | |
2766 | _dxfBuildSimplifiedSurface |
2767 | |
2768 +--------------------------------------------------------------------------+*/
2769
_dxfBuildSimplifiedSurface(SimpData * simp_data,int * new_nV,float ** new_v,int * new_nT,int ** new_t,float ** new_face_areas,float ** face_normals,float ** new_positional_error,float ** new_vertex_data)2770 int _dxfBuildSimplifiedSurface(SimpData *simp_data, int *new_nV, float **new_v,
2771 int *new_nT, int **new_t,
2772 float **new_face_areas, float **face_normals,
2773 float **new_positional_error, float **new_vertex_data)
2774
2775 /*...........................................................................*/
2776 /*
2777
2778 Use the information contained in the vertex and triangle fathers
2779 to construct a simplified surface:
2780
2781 only the vertices and triangles that are their own "father" are carried to
2782 the simplified surface. The vertex_parent and face_parent array are reshuffled
2783 to be used later for the purpose of interpolating data (positions-dependent
2784 or connections-dependent) from the original surface to the simplified surface)
2785
2786 I use indifferently the word "parent" or "father" in this
2787 implementation. I think "parent" is better than "father", which
2788 refers to the most "direct" parent. But the father of the father
2789 is a "parent", and it sounds better to me than a
2790 "father". However, there was much existing code using the word
2791 "father" so I kept it.
2792
2793
2794 */
2795 /*...........................................................................*/
2796
2797 {
2798
2799 /* count the new number of vertices */
2800
2801 int
2802 *vertex_parents = simp_data->vertex_father,
2803 *face_parents = simp_data->triangle_father;
2804
2805 register int i,j;
2806
2807 float
2808 *v = NULL;
2809
2810 int
2811 *t = NULL,
2812 *new_v_number = NULL,
2813 *new_t_number = NULL;
2814
2815
2816
2817 /* I) build the new vertices and copy the vertex errors and data */
2818
2819
2820
2821
2822 new_v_number = (int *) DXAllocateZero(simp_data->nV * sizeof (int));
2823
2824 if (!new_v_number) goto error;
2825
2826 *new_nV = 0;
2827
2828 for (i=0; i< simp_data->nV; i++) if (i == FATHER(i,vertex_parents))
2829
2830 {
2831 new_v_number[i] = *new_nV;
2832 *new_nV += 1;
2833 }
2834
2835
2836 /* allocate the new vertex array */
2837
2838 v = (float *) DXAllocateZero (*new_nV * sizeof (Vertex));
2839
2840 if (!v) goto error;
2841
2842 *new_positional_error = (float *) DXAllocateZero (*new_nV * sizeof (float));
2843
2844 if (!(*new_positional_error)) goto error;
2845
2846
2847 if (simp_data->vx_data)
2848 {
2849 *new_vertex_data = (float *) DXAllocateZero (*new_nV * simp_data->data_dim * sizeof (float));
2850
2851 if (!(*new_vertex_data)) goto error;
2852 }
2853
2854
2855 /* copy the vertices and other vertex attributes and update the vertex_parents array */
2856
2857
2858 for(j=i=0;i<simp_data->nV;i++) if (i == vertex_parents[i])
2859 {
2860 memcpy(v + 3 * j,simp_data->vert[i],sizeof(Vertex));
2861
2862 (*new_positional_error)[j] = simp_data->err_volume[i];
2863
2864 if (simp_data->vx_data)
2865 memcpy(*new_vertex_data + j * simp_data->data_dim,
2866 simp_data->vx_data + i * simp_data->data_dim, simp_data->data_dim * sizeof (float));
2867
2868 j++;
2869 }
2870
2871 for(i=0;i<simp_data->nV;i++) vertex_parents[i] = new_v_number[vertex_parents[i]];
2872
2873 DXFree((Pointer) new_v_number); new_v_number = NULL;
2874
2875
2876 /* copy the new vertices into *new_v and reallocate *new_v */
2877
2878
2879 DXFree((Pointer) (*new_v));
2880
2881 *new_v = (float *) DXAllocateZero(*new_nV * sizeof (Vertex));
2882
2883 if (!(*new_v)) goto error;
2884
2885 memcpy(*new_v, v, *new_nV * sizeof (Vertex));
2886
2887 DXFree((Pointer) v); v = NULL;
2888
2889
2890
2891
2892 /* II) build the new triangles and copy triangle normals and areas */
2893
2894
2895
2896
2897 new_t_number = (int *) DXAllocateZero(simp_data->nT * sizeof (int));
2898
2899 if (!new_t_number) goto error;
2900
2901 *new_nT = 0;
2902
2903 for (i=0; i< simp_data->nT; i++) if (i == FATHER(i,face_parents))
2904
2905 {
2906 new_t_number[i] = *new_nT;
2907 *new_nT += 1;
2908 }
2909
2910 t = (int *) DXAllocateZero (*new_nT * sizeof (Face));
2911
2912 if (!t) goto error;
2913
2914 *new_face_areas = (float *) DXAllocateZero (*new_nT * sizeof (float));
2915
2916 if (!(*new_face_areas)) goto error;
2917
2918 *face_normals = (float *) DXAllocateZero (*new_nT * sizeof (Vertex));
2919
2920 if (!(*face_normals)) goto error;
2921
2922 /* copy the triangles and other triangle attributes and update the face_parents array */
2923
2924 for(j=i=0;i<simp_data->nT;i++) if (i == face_parents[i])
2925 {
2926 register u_char k;
2927
2928 for (k=0;k<3;k++) t[j*3+k] = vertex_parents[ (*new_t)[i * 3 + k]];
2929
2930 (*new_face_areas)[j] = simp_data->area[i];
2931
2932 memcpy(*face_normals + 3 * j, simp_data->normal+i, sizeof(Vertex));
2933
2934 j++;
2935 }
2936
2937
2938 for(i=0;i<simp_data->nT;i++) face_parents[i] = new_t_number[face_parents[i]];
2939
2940
2941 DXFree((Pointer) new_t_number); new_t_number = NULL;
2942
2943
2944 /* copy the new triangles into *new_t and reallocate *new_t */
2945
2946 DXFree((Pointer) (*new_t));
2947
2948 *new_t = (int *) DXAllocateZero(*new_nT * sizeof (Face));
2949
2950 if (!(*new_t)) goto error;
2951
2952 memcpy(*new_t, t, *new_nT * sizeof (Face));
2953
2954 DXFree((Pointer) t); t = NULL;
2955
2956
2957 return 1;
2958
2959
2960 error:
2961
2962 if (new_t_number) DXFree((Pointer) new_t_number);
2963 if (new_v_number) DXFree((Pointer) new_v_number);
2964
2965 if (t) DXFree((Pointer) t);
2966 if (v) DXFree((Pointer) v);
2967
2968
2969 return 0;
2970 }
2971
2972 /*+--------------------------------------------------------------------------+
2973 | |
2974 | _dxfMarkEdgesAdjacent2Boundary |
2975 | |
2976 +--------------------------------------------------------------------------+*/
2977
2978
2979 /*
2980 As the name indicates, I am just using two different flags for edges adjacent to
2981 the surface boundary.
2982
2983 An edge that is a boundary edge, which is the same as an edge incident to two
2984 boundary vertex, is said to be a boundary edge of the second kind,
2985 with a flag value of -2
2986
2987 Otherwise, an edge incident to a single boundary vertex is said to
2988 be a boundary edge of the first kind, with a flag value of -1.
2989
2990 This routine is executed *only* if the user does not wish to simplify the boundary.
2991
2992 Once an edge is flagged with a negative flag value, it is impossible to add it
2993 later to the edge_heap.
2994
2995 */
2996
_dxfMarkEdgesAdjacent2Boundary(SimpData * simp_data)2997 int _dxfMarkEdgesAdjacent2Boundary(SimpData *simp_data)
2998 {
2999
3000 int
3001 e = 0,
3002 nE = simp_data->nE,
3003 edges_adjacent_2_boundary;
3004
3005 u_char
3006 *boundary_vertex = simp_data->boundary_vert,
3007 boundary1,
3008 boundary2;
3009
3010 EdgeS *current_edge = simp_data->edge_array;
3011
3012 /* for each edge of the surface, determine whether it is a boundary
3013 edge of the first or second type. If yes, label it accordingly
3014 using the edge2index array */
3015
3016 for (edges_adjacent_2_boundary = nE,
3017 e=0; e<nE; e++, current_edge++)
3018 {
3019
3020 /* if the edge label is 1, then it is a boundary edge of the second kind */
3021
3022 if (SurfEdgeLabel(current_edge) == 1) simp_data->edge2index[e] = -2;
3023
3024 else
3025 {
3026
3027 int
3028 *v = SurfEdgeGetVertices(current_edge);
3029
3030 boundary1 = boundary_vertex[v[0]];
3031 boundary2 = boundary_vertex[v[1]];
3032
3033 if (boundary1 || boundary2) simp_data->edge2index[e] = -1;
3034
3035 else edges_adjacent_2_boundary--;
3036 }
3037
3038 }
3039
3040
3041 return edges_adjacent_2_boundary;
3042 }
3043
3044 /*+--------------------------------------------------------------------------+
3045 | |
3046 | _dxfBuildEdgeHeap |
3047 | |
3048 +--------------------------------------------------------------------------+*/
3049
3050 /* add to the Heap edges all the edges that can be added */
3051
_dxfBuildEdgeHeap(SimpData * simp_data)3052 int _dxfBuildEdgeHeap(SimpData *simp_data)
3053 {
3054
3055 int
3056 can_add_edge = 1,
3057 n_edges_added = 0,
3058 nE = simp_data->nE;
3059
3060 while (can_add_edge && simp_data->last_edge_added2heap < nE)
3061 {
3062 if ( simp_data->edge2index[simp_data->last_edge_added2heap] ==
3063 simp_data->heap_initial_index )
3064
3065 /* *new edges* are entered at this stage. if an edge has
3066 been visited before, it can only be reentered if a
3067 neighbor was simplified. _dxfReinstateEdgesInHeap does this
3068 operation */
3069
3070 n_edges_added +=
3071 (can_add_edge =
3072 simp_data->add_edge(simp_data, simp_data->last_edge_added2heap));
3073
3074
3075 simp_data->last_edge_added2heap ++;
3076
3077 } /* end while */
3078
3079 return (n_edges_added);
3080 }
3081
3082 /*+--------------------------------------------------------------------------+
3083 | |
3084 | _dxfRotateParentTriangle |
3085 | |
3086 +--------------------------------------------------------------------------+*/
3087
3088
3089 /* performs a circular permutation (rotation) of the
3090 vertex indices of a triangle such that the parent of the
3091 first index is the same as "parent_vertex" */
3092
_dxfRotateParentTriangle(int parent_vertex,int * vertex_fathers,int * tri_v,int * tri_v_parents,int * tri_v_rotated)3093 int _dxfRotateParentTriangle(int parent_vertex, int *vertex_fathers,
3094 int *tri_v, int *tri_v_parents, int *tri_v_rotated)
3095 {
3096
3097 /* the triangle must be a parent triangle itself (otherwise
3098 all three triangle vertices are not necessarily different */
3099
3100 tri_v_parents[0] = FATHER(tri_v[0], vertex_fathers);
3101 tri_v_parents[1] = FATHER(tri_v[1], vertex_fathers);
3102 tri_v_parents[2] = FATHER(tri_v[2], vertex_fathers);
3103
3104 if (tri_v_parents[0] == parent_vertex)
3105 {
3106 tri_v_rotated[0] = tri_v[0];
3107 tri_v_rotated[1] = tri_v[1];
3108 tri_v_rotated[2] = tri_v[2];
3109 }
3110 else if (tri_v_parents[1] == parent_vertex)
3111 {
3112 /* rotate the vertex parents */
3113
3114 tri_v_rotated[0] = tri_v_parents[0]; /* use tri_v_rotated[0] as a temporary buffer */
3115
3116 tri_v_parents[0] = tri_v_parents[1];
3117 tri_v_parents[1] = tri_v_parents[2];
3118 tri_v_parents[2] = tri_v_rotated[0];
3119
3120 /* rotate the vertices */
3121 tri_v_rotated[0] = tri_v[1];
3122 tri_v_rotated[1] = tri_v[2];
3123 tri_v_rotated[2] = tri_v[0];
3124 }
3125 else
3126 {
3127 /* rotate the vertex parents */
3128 tri_v_rotated[0] = tri_v_parents[0]; /* use tri_v_rotated[0] as a temporary buffer */
3129
3130 tri_v_parents[0] = tri_v_parents[2];
3131 tri_v_parents[2] = tri_v_parents[1];
3132 tri_v_parents[1] = tri_v_rotated[0];
3133
3134 /* rotate the vertices */
3135 tri_v_rotated[0] = tri_v[2];
3136 tri_v_rotated[1] = tri_v[0];
3137 tri_v_rotated[2] = tri_v[1];
3138 }
3139
3140 return 1;
3141 }
3142
3143 /*+--------------------------------------------------------------------------+
3144 | |
3145 | _dxfRotateParentEdge |
3146 | |
3147 +--------------------------------------------------------------------------+*/
3148
_dxfRotateParentEdge(int parent_vertex,int * vertex_fathers,int * edge_v,int * edge_t,int * edge_v_parents,int * edge_v_rotated,int * edge_t_rotated)3149 int _dxfRotateParentEdge(int parent_vertex, int *vertex_fathers,
3150 int *edge_v, int *edge_t, int *edge_v_parents,
3151 int *edge_v_rotated, int *edge_t_rotated)
3152 {
3153
3154
3155 /* rotate a parent edge with respect to a parent vertex, so that
3156 the child of the parent vertex will be first (index [0])
3157 among the two edge vertices
3158 ( currently, by default the first vertex is the vertex with the
3159 highest number )
3160 */
3161
3162 edge_v_parents[0] = FATHER(edge_v[0], vertex_fathers);
3163 edge_v_parents[1] = FATHER(edge_v[1], vertex_fathers);
3164
3165 /* if the vertices are in the right order, just copy the vertices */
3166
3167 if (edge_v_parents[0] == parent_vertex)
3168 {
3169 edge_v_rotated[0] = edge_v[0];
3170 edge_v_rotated[1] = edge_v[1];
3171 edge_t_rotated[0] = edge_t[0];
3172 edge_t_rotated[1] = edge_t[1];
3173 }
3174 else
3175 {
3176 edge_v_rotated[0] = edge_v_parents[0]; /* use edge_v_rotated[0] as a temporary buffer */
3177
3178 edge_v_parents[0] = edge_v_parents[1];
3179 edge_v_parents[1] = edge_v_rotated[0];
3180
3181 edge_v_rotated[0] = edge_v[1];
3182 edge_v_rotated[1] = edge_v[0];
3183 edge_t_rotated[0] = edge_t[1];
3184 edge_t_rotated[1] = edge_t[0];
3185
3186 }
3187
3188 return 1;
3189 }
3190
3191 /*+--------------------------------------------------------------------------+
3192 | |
3193 | _dxfDirectedParentVertexStar |
3194 | |
3195 +--------------------------------------------------------------------------+*/
3196
3197 /*
3198 compute a vertex star and store the triangles in either
3199 CLOCKWISE or COUNTER-CLOCKWISE direction
3200 */
3201
_dxfDirectedParentVertexStar(int parent_vertex,int first_triangle,int valence,int * star,SimpData * simp_data,int direction)3202 int _dxfDirectedParentVertexStar(int parent_vertex, int first_triangle,
3203 int valence, int *star, SimpData *simp_data, int direction)
3204 {
3205
3206
3207 /* labeling system for directed triangle stars
3208
3209 v*1\ /v*0
3210 e*1\ *1/e*0
3211 *val0-1 \ / *0
3212 ___________\/___________ we omit on purpose the first vertex, because in the case
3213 v*val0-1 v0 of a boundary edge, we know it already
3214
3215
3216 */
3217
3218 int
3219 *vstar = star + valence,
3220 *estar = vstar + valence,
3221 /* rotated triangle vertices: */
3222 tri_v_rotated[3], tri_v_parents[3],
3223 /* rotated edge vertices and adjacent triangles: */
3224 edge_v_rotated[2], edge_v_parents[2], edge_t_rotated[2],
3225 star_triangle = first_triangle,
3226 i;
3227
3228
3229 for (i = 0; i < valence; i++)
3230 {
3231 star[i] = star_triangle;
3232
3233 _dxfRotateParentTriangle(parent_vertex, simp_data->vertex_father,
3234 simp_data->triangles + 3 * star_triangle,
3235 tri_v_parents, tri_v_rotated);
3236
3237 /* find previous or next edges and get their parents */
3238
3239 estar[i] = _dxfNextParentEdge(tri_v_rotated, direction);
3240
3241 _dxfRotateParentEdge(parent_vertex, simp_data->vertex_father,
3242 (int *) SurfEdgeGetVertices (simp_data->edge_array + estar[i]),
3243 (int *) SurfEdgeGetTriangles(simp_data->edge_array + estar[i]),
3244 edge_v_parents, edge_v_rotated, edge_t_rotated);
3245
3246 vstar[i] = edge_v_parents[1];
3247 /* this vertex number should be the same as tri_v_parents[direction] */
3248
3249 if (i < valence-1)
3250 /* otherwise, the next triangle would be -1 and finding its father
3251 would produce a segmentation fault or another bad thing */
3252 star_triangle = _dxfNextParentTriangle(edge_t_rotated, direction);
3253 }
3254
3255
3256 return 1;
3257 }
3258
3259 /*+--------------------------------------------------------------------------+
3260 | |
3261 | _dxfManifoldLinkTest |
3262 | |
3263 +--------------------------------------------------------------------------+*/
3264
3265 /*
3266 detect duplicate vertices in the link of an edge.
3267 In case duplicate vertices are detected, collapsing an edge
3268 would generate a singular vertex.
3269 */
3270
_dxfManifoldLinkTest(int * link0,int * link1,int val0,int val1)3271 int _dxfManifoldLinkTest(int *link0, int *link1, int val0, int val1)
3272 {
3273 /* among the vertices of the link of v0, link0 and of the link of v1, link1,
3274 is there any duplicate vertex ? */
3275
3276 int
3277 link_is_manifold = 1,
3278 val = val0 + val1;
3279
3280 /* since
3281 val0 + val1 <= SIMPLIFY_VALENCE_MAX4,
3282 I can use a static array here. */
3283
3284 static int
3285 link[SIMPLIFY_VALENCE_MAX4];
3286
3287 {
3288 register int i;
3289
3290 /* gather vertices of link0 and link1 in a single array, sort
3291 the array, and determine whether there are any duplicate
3292 vertices */
3293
3294 memcpy(link, link0, val0 * sizeof (int));
3295 memcpy(link + val0, link1, val1 * sizeof (int));
3296
3297 qsort(link, val, sizeof(int), (int (*)(const void*, const void*)) _dxfCmpIntSmallFirst);
3298
3299 i = 0;
3300
3301 while ((i < val -1 ) && (link_is_manifold))
3302 {
3303 if (link[i+1] == link[i]) /* if there are duplicate vertices,
3304 they must be contiguous in the
3305 sorted array */
3306
3307 link_is_manifold = 0;
3308 i++;
3309 }
3310
3311 }
3312
3313 return link_is_manifold;
3314
3315 }
3316
3317 /*+--------------------------------------------------------------------------+
3318 | |
3319 | _dxfCmpIntSmallFirst |
3320 | |
3321 +--------------------------------------------------------------------------+*/
3322
3323 /*
3324 compare integers. Using this function in qsort will
3325 result in listing the smaller integers first
3326 */
3327
_dxfCmpIntSmallFirst(char * s1,char * s2)3328 int _dxfCmpIntSmallFirst(char *s1, char *s2)
3329 {
3330 int
3331 *f1 = (int *)s1,
3332 *f2 = (int *)s2;
3333
3334 if (*f2 > *f1)
3335 return(-1);
3336
3337 else if (*f2 < *f1)
3338 return(1);
3339
3340 else
3341 return(0);
3342
3343 }
3344 /*+--------------------------------------------------------------------------+
3345 | |
3346 | _dxfMakeEdgeBarycenter |
3347 | |
3348 +--------------------------------------------------------------------------+*/
3349
3350 /*
3351 compute the barycenter of two vertices
3352 */
3353
_dxfMakeEdgeBarycenter(Vertex v0,Vertex v1,float alpha0,Vertex v)3354 void _dxfMakeEdgeBarycenter(Vertex v0, Vertex v1, float alpha0, Vertex v)
3355 {
3356
3357 register u_char j;
3358
3359 float alpha1 = 1. - alpha0;
3360
3361 for(j=0; j<3; j++)
3362
3363 v[j] = alpha0 * v0[j] + alpha1 * v1[j];
3364
3365 return;
3366 }
3367
3368 /*+--------------------------------------------------------------------------+
3369 | |
3370 | _dxfClosestPointOnEdge |
3371 | |
3372 +--------------------------------------------------------------------------+*/
3373
3374 /*
3375 compute the closest point XP to a point X on an edge (A,B)
3376
3377 return the barycentric cordinnates of X with respect to A and B
3378
3379 */
3380
3381
_dxfClosestPointOnEdge(Vertex X,Vertex XP,Vertex A,Vertex B,float * t)3382 float _dxfClosestPointOnEdge(Vertex X, Vertex XP, Vertex A, Vertex B,
3383 float *t)/* barycentric coordinates: t, 1-t */
3384 {
3385 float d;
3386 Vertex BA,BX;
3387 float BA_BX,BA_2;
3388
3389 MakeVec(B,A,BA);
3390 MakeVec(B,X,BX);
3391 EdgeBarycentricCoordinates(BA,BX,BA_BX,BA_2);
3392
3393 if (BA_BX <= 0.)
3394 {
3395 *t = 0.;
3396 d = NORM3(BX);
3397 VecCopy(B,XP);
3398 }
3399 else if (BA_2 <= BA_BX)
3400 {
3401 *t = 1.;
3402 d = DIST3(A,X);
3403 VecCopy(A,XP);
3404 }
3405 else
3406 {
3407 /* the following is not robust numerically */
3408 /*
3409 d= SCALPROD3(BX,BX) - BA_BX*BA_BX / BA_2;
3410 if (d<=0.0) d = 0.0;
3411 else d =(float) sqrt ((double) d);
3412 */
3413
3414 /* it is more robust to compute d= NORM3(Vecprod(BA,BC))/BA
3415 or to just use the distance from X to XP */
3416
3417 *t = BA_BX/BA_2;
3418
3419 _dxfMakeEdgeBarycenter(A,B,*t,XP);
3420
3421 d = DIST3(X,XP);
3422
3423 }
3424
3425
3426 return(d);
3427 }
3428
3429 /*+--------------------------------------------------------------------------+
3430 | |
3431 | _dxfNormalize3Vector |
3432 | |
3433 +--------------------------------------------------------------------------+*/
3434
_dxfNormalize3Vector(Vertex v)3435 float _dxfNormalize3Vector(Vertex v)
3436 {
3437 register u_char j;
3438 register double norm;
3439
3440 for (norm=0.0,j=0;j<3;j++) norm+=v[j]*v[j];
3441
3442
3443 if (norm>0.0)
3444 {
3445 norm=sqrt(norm);
3446 for (j=0;j<3;j++) v[j] /= norm;
3447 }
3448
3449 return (float)norm ;
3450 }
3451
3452 /*+--------------------------------------------------------------------------+
3453 | |
3454 | _dxfSolveQuadraticEquation |
3455 | |
3456 +--------------------------------------------------------------------------+*/
3457
_dxfSolveQuadraticEquation(double A,double B,double C,double * sol1,double * sol2,double eps)3458 int _dxfSolveQuadraticEquation(double A, double B, double C,
3459 double *sol1, double *sol2, double eps)
3460 {
3461 /*
3462
3463 returns the number of real solutions of a second degree equation as well as the solutions
3464
3465 everything below eps is supposed to be 0
3466
3467 */
3468
3469 int num_sol;
3470
3471 A *= 2.;
3472
3473 if ( fabs( A ) < eps)
3474 {
3475 /* B x + C = 0 */
3476
3477 if ( fabs( B ) < eps)
3478 num_sol = 0;
3479
3480 else
3481 {
3482 num_sol = 1;
3483 *sol1 = - C / B;
3484 }
3485 }
3486
3487 else
3488 {
3489 double delta = B * B - 2. * A * C;
3490
3491 if ( delta < - eps )
3492 num_sol = 0;
3493
3494 else
3495 {
3496 /* the discriminant is zero or positive */
3497
3498 if (delta < eps)
3499 {
3500 num_sol = 1;
3501
3502 *sol1 = -B / A;
3503 }
3504 else
3505 {
3506 delta = sqrt (delta);
3507
3508 num_sol = 2;
3509
3510 *sol1 = (-B - delta ) /A;
3511
3512 *sol2 = (-B + delta ) /A;
3513 }
3514 }
3515 }
3516
3517 return num_sol;
3518
3519 }
3520 /*+--------------------------------------------------------------------------+
3521 | |
3522 | _dxfSimplifiedBoundaryVertexLocation |
3523 | |
3524 +--------------------------------------------------------------------------+*/
3525
3526 /*
3527 determine the location of a potential simplified boundary vertex.
3528 applies only when a boundary edge is collapsed, which implies that
3529 the two endpoints of the edge are boundary vertices
3530 */
3531
3532
_dxfSimplifiedBoundaryVertexLocation(SimpData * simp_data,int v0,int v1,int * vstar0,int * vstar1,int val0,int val1,int method)3533 int _dxfSimplifiedBoundaryVertexLocation(SimpData *simp_data, int v0, int v1,
3534 int *vstar0, int *vstar1,
3535 int val0, int val1, int method)
3536 {
3537 int success = 1;
3538
3539 if (method == 0)
3540
3541 memcpy(simp_data->simplified_vertex, simp_data->vert[v0], sizeof (Vertex));
3542
3543 else if (method == 1)
3544 {
3545
3546 /* place the simplified vertex so that the length of the boundary
3547 will be preserved. The vertex must be situated on an ellipse
3548
3549 ___________\/___________\/__________
3550 v*val0-1 v0 v1 v*val1-1
3551 == vl == vr
3552
3553 Firstly, we retrieve the vertices vl, v0, v1, vr from the stars.
3554
3555 we compute the length of the polygonal arc vl--v0--v1--vr, and
3556 call 2a this length
3557
3558 The focal points of the ellipsoid are vr and vl, so we note 2c the distance
3559 vr-vl.
3560
3561 */
3562
3563 float
3564 a, a2, c, eps1 = (float)1e-6,
3565 *ve = (float *) simp_data->simplified_vertex;
3566
3567 Vertex
3568 v_zero, v_one, vr, vl, vm;
3569
3570 memcpy(v_zero, simp_data->vert[v0], sizeof (Vertex));
3571 memcpy(v_one, simp_data->vert[v1], sizeof (Vertex));
3572
3573 memcpy(vl, simp_data->vert[vstar0[val0-1]], sizeof (Vertex));
3574 memcpy(vr, simp_data->vert[vstar1[val1-1]], sizeof (Vertex));
3575
3576 /* compute the parameters of the ellipsoid */
3577
3578 a =
3579 (
3580 DIST3(simp_data->vert[v0], simp_data->vert[v1]) +
3581 DIST3(simp_data->vert[vstar0[val0-1]], simp_data->vert[v0]) +
3582 DIST3(simp_data->vert[vstar1[val1-1]], simp_data->vert[v1])
3583 )
3584 / 2.;
3585
3586 a2 = a * a;
3587
3588 c = DIST3(vl,vr)/2.;
3589
3590 /* Compute the middle of v_zero and v_one */
3591
3592 vm[0] = (v_zero[0] + v_one[0]) / 2.;
3593 vm[1] = (v_zero[1] + v_one[1]) / 2.;
3594 vm[2] = (v_zero[2] + v_one[2]) / 2.;
3595
3596 /* Then, our first next test is to determine whether a and c are
3597 sufficiently different, otherwise, taking the middle of v_zero and
3598 v_one will certainly work since we will have the same cord length
3599 along the boundary */
3600
3601 if (a == 0.0 || /* if a == 0, the length to preserve is zero, then the segment middle is fine! */
3602 1. - c/a < eps1)
3603
3604 /* limit of the floating point precision */
3605
3606 {
3607 memcpy(ve, vm, sizeof(Vertex));
3608 }
3609
3610 else /* if (a > c ) */
3611 {
3612 float
3613 t, x, y,
3614 ye;
3615
3616 Vertex
3617 vp;
3618
3619 /* then we consider the middle of v0 and v1, vm,
3620 we project vm on the segment vr vl and obtain vp
3621
3622 finally, we project vp on the ellipsoid defined by vr, vl, a
3623 and c along the direction vp--->vm
3624
3625
3626 ve
3627
3628
3629 vm
3630 |y
3631 vl________________0____|_____________vr
3632 -c x c
3633 vp
3634
3635
3636 by projecting vm onto the segment (vl,vc), we obtain
3637 vp, and the coordinates of vm in the coordinate system formed by
3638 0(the middle of vr and vl) vr and vm.
3639
3640 x and y are the coordinates of vm in this coordinate system
3641 We keep x unchanged, and compute ye such that
3642
3643 x^2/a^2 + ye^2/b^2 = 1, with b^2 = a^2-c^2
3644
3645 x^2
3646 hence ye^2 = ( 1 - --- ) ( a^2 - c^2 )
3647 a^2
3648
3649 We take for ye the positive square root of ye^2
3650
3651 and finally ve is obtained by
3652
3653 ve-vp = ye/y (vm-vp)
3654
3655 ve = vp + ye/y (vm-vp)
3656
3657
3658 5/7/97 bug of Floating Point exception noticed on DEC
3659 A special case is encountered if c == 0,
3660 then we also have x = 0 and the formula is simply ye=a
3661 this case is handled with the general case without particular
3662 attention, except that we can't divide by c and
3663 also we should not execute the part of the code relative to t == 0 or t == 1
3664 */
3665
3666
3667 /* compute x and y of vm with respect to vr and vl */
3668
3669 y = _dxfClosestPointOnEdge(vm, vp, vl, vr, &t);
3670
3671
3672 /* what do we do if y == 0, in which case vm is accidentally
3673 exactly on the line segment vr-vl?
3674
3675 In that case, we still have a non-degenerate ellipse since we
3676 ruled out the fact that a and c are the same. Then either
3677 v_zero and v_one cannot be on the line segment vr-vl.
3678
3679 So we project instead v0 on the segment vr vl
3680
3681 v0
3682 ___|________________
3683 vl vp0 vp vr
3684
3685 ----->
3686 and we use the direction vp0 v0 for projecting the point vp on the ellipse
3687
3688 */
3689
3690
3691 if (y <= eps1 * c) /* it can be that c == 0 and y == 0 */
3692 {
3693 /* limit of floating point precision */
3694 /* In that case, either v_zero or v_one cannot be on the line vl, vr,
3695 so we replace vm with v_zero */
3696
3697 memcpy(vm, v_zero, sizeof(Vertex));
3698
3699 y = _dxfClosestPointOnEdge(vm, vp, vl, vr, &t);
3700
3701 if (y <= eps1 *c)
3702 {
3703 memcpy(vm, v_one, sizeof(Vertex));
3704
3705 y = _dxfClosestPointOnEdge(vm, vp, vl, vr, &t);
3706 }
3707 }
3708
3709 /* what would it mean for y to be equal to zero after these
3710 various projections: that all points v_zero, v_one, vm are on vl_vr,
3711 meaning that a = c. Just in case, we use an if (y!=0.0) in the end */
3712
3713 /* in case vm does not project orthogonally on the segment (vl,vr)
3714 the routine ClosestPointOnEdge projects vm either on vl and assigns
3715 1 to t or on vr and assigns 0 to t
3716 we need to determine the angle that vm,vp makes with (vl,vr)
3717 */
3718
3719 if (c > 0.0 /* c == 0 is handled as the general case */
3720 && (t == 0. || t == 1.))
3721 {
3722 Vertex u,v;
3723 double
3724 cos_uv, cos_uv_2, sin_uv_2,
3725 A, B, C, sol1, sol2, eps2 = 1e-12;
3726
3727 MakeVec(vl,vr,u); /* if the segment (vr,vl) is of zero length, which is
3728 possible in nasty cases, u will be a vector of zero
3729 lenght, and this invalidates the following computations
3730 we need to find another solution */
3731 _dxfNormalize3Vector(u);
3732
3733 MakeVec(vp,vm,v);
3734 _dxfNormalize3Vector(v);
3735
3736 cos_uv = u[0] * v[0] + u[1]*v[1] + u[2]*v[2]; /* FScalProd(u,v); */
3737
3738 if (cos_uv < 0.)
3739 cos_uv = -cos_uv;
3740
3741 cos_uv_2 = cos_uv * cos_uv;
3742
3743 sin_uv_2 = 1. - cos_uv_2;
3744
3745 /* ye is the solution of the following equation of the second
3746 degree:
3747
3748 x^2/a^2 + y^2/(a^2-c^2) = 1 <=>
3749
3750 (ye cos_uv + c)^2/a^2 + (ye sin_uv)^2/(a^2-c^2) = 1 <=>
3751
3752 ye^2( cos_uv_2 + a^2/(a^2-c^2) sin_uv_2) +
3753 ye (2 c cos_uv) +
3754 c^2-a^2 = 0
3755
3756
3757 A ye^2 + cb ye + C = 0 */
3758
3759 A = cos_uv_2 + sin_uv_2 * a2 / (a+c)/(a-c);
3760 B = 2. * cos_uv * c;
3761 C = (c+a)*(c-a); /* c^2-a^2 */
3762
3763 _dxfSolveQuadraticEquation(A, B, C, &sol1, &sol2, eps2);
3764
3765 /* there is only one positive solution */
3766
3767 if (sol1 > 0.)
3768 ye = sol1;
3769 else if (sol2 > 0.)
3770 ye = sol2;
3771 else
3772 {
3773 ye = 0;
3774
3775 /* The quadratic equation has no positive root */
3776 }
3777 } /*endif vm project on an extremity */
3778
3779 else /* vm projects normally inside the segment */
3780 {
3781 double tmp;
3782
3783 /* determine x from t: */
3784
3785 /* vp = t * vl + (1-t) * vr
3786 t == 1 <==> x = -c
3787 t == 0 <==> x = c
3788
3789 x = c - 2 c t = c ( 1-2t) */
3790
3791 x = c * (1. - 2. * t);
3792
3793
3794 /* ye = sqrt ((double) (a2 - c * c)*( 1. - x*x/a2)); */
3795
3796 /* since a and c are likely to be very close, a more robust
3797 computation is to avoid computations such as a2 - c^2
3798 and compute sqrt(a2) is a is readily available */
3799
3800 tmp = (double) (a -c) * (a+c) * (a-x) * (a+x);
3801
3802 /* this is to avoid floating point exceptions when computing the square root */
3803
3804 if (tmp >0.0) ye = sqrt (tmp) / a;
3805
3806 else ye = 0.0;
3807
3808 }
3809
3810 if (y != 0.0)
3811 ye /= y;
3812
3813 /* do the following vector composition without
3814 a function call ve = vp + ye/y (vm-vp) */
3815
3816 ve[0] = vp[0] + ye * (vm[0] - vp[0]);
3817 ve[1] = vp[1] + ye * (vm[1] - vp[1]);
3818 ve[2] = vp[2] + ye * (vm[2] - vp[2]);
3819
3820 /* test whether the boundary length is preserved */
3821 {
3822
3823 float
3824 l_before = 2. * a,
3825 l_after = DIST3(vl,ve) + DIST3(ve,vr),
3826
3827 eps3 = simp_data->tolerance * eps1; /* this is a relative accuracy
3828 tolerance * floating point accuracy
3829 using eps1 here would not be general enough */
3830
3831 if (l_before -l_after > eps3 || l_before - l_after < -eps3)
3832
3833 {
3834 success = 0;
3835 }
3836
3837 }
3838
3839 } /* if (a>c) otherwise, vm, the middle is good for ve */
3840
3841 } /* if (method = 1, length preservation) otherwise */
3842
3843 return success;
3844 }
3845
3846 /*+--------------------------------------------------------------------------+
3847 | |
3848 | _dxfFastCompactness3DTriangle2 |
3849 | |
3850 +--------------------------------------------------------------------------+*/
3851
3852 /* compute the triangle compactness and area
3853 this is a fast method, in comparison with the method that
3854 orthogonalizes the triangle and computes the area afterwards
3855
3856 this fast method uses a cross product
3857 */
3858
_dxfFastCompactness3DTriangle2(Vertex * tri,float * triangle_area)3859 float _dxfFastCompactness3DTriangle2(Vertex *tri, float *triangle_area)
3860 {
3861
3862 /* compute the squared lenghts of the three edges: */
3863
3864 double sum_lengths_sq = 0., tmp;
3865
3866 {
3867
3868 register u_char i,j;
3869
3870 for(i=0;i<3;i++)
3871 {
3872 double len_edg_sq = 0;
3873
3874 for(j=0;j<3;j++)
3875 {
3876 tmp = tri[(i+1)%3][j] - tri[i][j];
3877 len_edg_sq += tmp * tmp;
3878 }
3879
3880 sum_lengths_sq += len_edg_sq;
3881
3882 }
3883 }
3884
3885
3886 {
3887 double
3888 m0 =
3889 tri[1][1] - tri[0][1],
3890 m1 =
3891 tri[2][2] - tri[1][2],
3892 m2 =
3893 tri[1][2] - tri[0][2],
3894 m3 =
3895 tri[2][1] - tri[1][1],
3896 m4 =
3897 tri[2][0] - tri[1][0],
3898 m5 =
3899 tri[1][0] - tri[0][0],
3900 n0 =
3901 m0 * m1 -
3902 m2 * m3,
3903
3904 n1 =
3905 m2 * m4 -
3906 m5 * m1,
3907
3908 n2 =
3909 m5 * m3 -
3910 m0 * m4;
3911
3912 tmp = n0 * n0 + n1 * n1 + n2 * n2;
3913
3914 if (tmp > 0.0) *triangle_area = sqrt (tmp) / 2.;
3915
3916 else *triangle_area = 0;
3917
3918
3919 }
3920
3921
3922 return (float) (FOUR_SQRT_3 * *triangle_area / sum_lengths_sq);
3923
3924 }
3925
3926
3927 /*+--------------------------------------------------------------------------+
3928 | |
3929 | _dxfBoundaryCollapse2KndGeometricallyOk |
3930 | |
3931 +--------------------------------------------------------------------------+*/
3932
_dxfBoundaryCollapse2KndGeometricallyOk(SimpData * simp_data,int v0,int v1,int * star0,int * vstar0,int val0,Vertex * s_normal0,float * s_area0,int direction,float min_scalprod_nor,float compactness_ratio)3933 int _dxfBoundaryCollapse2KndGeometricallyOk(SimpData *simp_data, int v0, int v1,
3934
3935 /* this procedure must be run twice, once for the
3936 star of the vertex v0 and once for the star of the vertex v1
3937 Here, v1 rather stands for the first vertex of the star
3938 */
3939
3940 int *star0, int *vstar0, int val0, Vertex *s_normal0, float *s_area0,
3941 int direction, float min_scalprod_nor, float compactness_ratio)
3942
3943 {
3944 int collapsible = 1, i = 0;
3945
3946 Vertex faceVert[3];
3947
3948 float
3949 *s_comp0 = s_area0 + val0,
3950 min_compactness_before,
3951 min_compactness_after;
3952
3953 /* labeling system for directed triangle stars
3954
3955 v*1\ /v*0
3956 e*1\ *1/e*0
3957 *val0-1 \ / *0
3958 ___________\/___________ we omit on purpose the first vertex, because in the case
3959 v*val0-1 v0 v1 of a boundary edge, we know it already
3960
3961
3962 */
3963
3964 /* retrieve the compactness of the first triangle */
3965
3966 min_compactness_before = simp_data->compactness[star0[0]];
3967
3968 min_compactness_after = 1; /* initialization of the compactness
3969 after simplification */
3970
3971 while (collapsible && i < val0-1)
3972 {
3973
3974 /* retrieve the face compactness before simplification */
3975
3976 min_compactness_before = MIN ( simp_data->compactness[star0[i+1]], min_compactness_before);
3977
3978 /* recompute the compactness after changing the first vertex */
3979
3980 memcpy(faceVert[0], simp_data->simplified_vertex, sizeof (Vertex));
3981
3982 if (direction == DIRECTION_CLOCKWISE)
3983 {
3984 /* reverse the order of triangle vertices, so that the normal
3985 that we will compute will respect the surface orientation */
3986 memcpy(faceVert[2], simp_data->vert[vstar0[i]], sizeof (Vertex));
3987 memcpy(faceVert[1], simp_data->vert[vstar0[i+1]], sizeof (Vertex));
3988 }
3989
3990 else
3991 {
3992 memcpy(faceVert[1], simp_data->vert[vstar0[i]], sizeof (Vertex));
3993 memcpy(faceVert[2], simp_data->vert[vstar0[i+1]], sizeof (Vertex));
3994 }
3995
3996 s_comp0[i+1] = _dxfFastCompactness3DTriangle2(faceVert, s_area0+i+1);
3997
3998 min_compactness_after = MIN ( s_comp0[i+1], min_compactness_after);
3999
4000 /* recompute the face normal assuming that the vertex v0 has
4001 moved and save this normal for future use*/
4002
4003 _dxfTriangleNormalQR2(faceVert, s_normal0[i+1]);
4004
4005
4006 /* check whether the triangle normal orientation is consistent */
4007
4008 collapsible = ((SCALPROD3(s_normal0[i+1], simp_data->normal[star0[i+1]])) > min_scalprod_nor);
4009
4010 i++;
4011 }
4012
4013 if (collapsible)
4014
4015 /* check whether the minimum triangle compactness is too much degraded */
4016
4017 collapsible = min_compactness_after > (compactness_ratio * min_compactness_before);
4018
4019 return collapsible;
4020 }
4021
4022 /*+--------------------------------------------------------------------------+
4023 | |
4024 | _dxfErrorTestDataOnSurfaceUpdate |
4025 | |
4026 +--------------------------------------------------------------------------+*/
4027
4028 /*
4029 given a new vertex of the mapping between simplified and original surfaces,
4030 update the potential data error value at the potential simplified vertex location
4031
4032 */
4033
4034
4035
_dxfErrorTestDataOnSurfaceUpdate(SimpData * simp_data,int new_v,int old_v,int new_v2,int new_v3,float alpha_sv,float alpha_n_v2,float alpha_n_v3,int old_v1,int old_v2,int old_v3,float alpha_o_v1,float alpha_o_v2,float alpha_o_v3)4036 int _dxfErrorTestDataOnSurfaceUpdate(SimpData *simp_data, int new_v, int old_v,
4037 int new_v2, int new_v3,
4038 float alpha_sv, float alpha_n_v2, float alpha_n_v3,
4039 int old_v1, int old_v2, int old_v3,
4040 float alpha_o_v1, float alpha_o_v2, float alpha_o_v3)
4041 {
4042
4043 int do_collapse = 1;
4044
4045 float
4046 o_alpha[3],
4047 n_alpha[3],
4048 *o_data[3],
4049 *n_data[3],
4050 o_data_error[3],
4051 data_difference,
4052 new_error_link,
4053 tmp_data_diff=0,
4054 new_error_sv,
4055 old_error;
4056
4057 register int j;
4058 register u_char k;
4059
4060 o_data[0] = simp_data->vx_data + old_v1 * simp_data->data_dim;
4061 o_data[1] = simp_data->vx_data + old_v2 * simp_data->data_dim;
4062 o_data[2] = simp_data->vx_data + old_v3 * simp_data->data_dim;
4063
4064 n_data[0] = simp_data->vx_data_potential_values;
4065 n_data[1] = simp_data->vx_data + new_v2 * simp_data->data_dim;
4066 n_data[2] = simp_data->vx_data + new_v3 * simp_data->data_dim;
4067
4068 o_alpha[0] = alpha_o_v1;
4069 o_alpha[1] = alpha_o_v2;
4070 o_alpha[2] = alpha_o_v3;
4071
4072 n_alpha[0] = alpha_sv;
4073 n_alpha[1] = alpha_n_v2;
4074 n_alpha[2] = alpha_n_v3;
4075
4076 o_data_error[0] = simp_data->vx_data_error[old_v1];
4077 if (old_v > 1)
4078 o_data_error[1] = simp_data->vx_data_error[old_v2];
4079 if (old_v > 2)
4080 o_data_error[2] = simp_data->vx_data_error[old_v3];
4081
4082
4083
4084 /* interpolate vertex data in the old and new surfaces
4085 and decide if they are sufficiently similar */
4086
4087 for (data_difference = 0., j=0; j<simp_data->data_dim; j++)
4088 {
4089 simp_data->vx_old_data[j] = 0.0;
4090
4091 for (k=0 ; k < old_v; k++) simp_data->vx_old_data[j] += o_alpha[k] * o_data[k][j];
4092
4093 simp_data->vx_new_data[j] = 0.0;
4094
4095 for (k=0 ; k < new_v; k++) simp_data->vx_new_data[j] += n_alpha[k] * n_data[k][j];
4096
4097 tmp_data_diff = simp_data->vx_new_data[j] - simp_data->vx_old_data[j];
4098
4099 data_difference += tmp_data_diff * tmp_data_diff;
4100
4101 }
4102
4103
4104 /* there are two cases, either the data dimension is one or it is larger than one
4105 in the first case, no square root is required for computing the discrepancy
4106 between data values, and since I am very stingy about expensive
4107 floating point computations, I am avoiding taking this square root if I can.
4108 */
4109
4110 if (simp_data->data_dim == 1) {data_difference = (tmp_data_diff > 0.0) ? tmp_data_diff : -tmp_data_diff;}
4111
4112 else
4113 if (data_difference > 0.0) data_difference = sqrt (data_difference);
4114
4115 old_error = 0.0;
4116
4117 for (k=0 ; k < old_v; k++) old_error += o_alpha[k] * o_data_error[k];
4118
4119 /* can the sum of the data difference and the old error at the
4120 current mesh point be explained by a new error at the simplified
4121 vertex, and if it is the case, what is this new error ? */
4122
4123 new_error_link = (new_v == 3) ?
4124 alpha_n_v2 * simp_data->vx_data_error[new_v2] + alpha_n_v3 * simp_data->vx_data_error[new_v3]:
4125 (new_v == 2) ? alpha_n_v2 * simp_data->vx_data_error[new_v2]:
4126 0.0;
4127
4128
4129 /* there is a special case if the simplified vertex is associated with an
4130 alpha that's zero: in this case, the error must be explained with the error at
4131 the link */
4132
4133 if (alpha_sv == 0. )
4134 {
4135 if (new_error_link < old_error + data_difference)
4136 do_collapse = 0;
4137
4138 /* but otherwise no potential error is computed */
4139
4140 }
4141 else /* alpha_sv * error_sv + new_error_link > old_error + data_difference */
4142 {
4143 new_error_sv = (old_error + data_difference - new_error_link) / alpha_sv;
4144
4145 if (new_error_sv > simp_data->vx_data_potential_err) simp_data->vx_data_potential_err = new_error_sv;
4146
4147 do_collapse = (simp_data->vx_data_potential_err <= simp_data->data_tolerance);
4148 }
4149
4150
4151 return do_collapse;
4152 }
4153
4154 /*+--------------------------------------------------------------------------+
4155 | |
4156 | _dxfSegmentIntersection |
4157 | |
4158 +--------------------------------------------------------------------------+*/
4159
_dxfSegmentIntersection(Vertex A,Vertex B,Vertex C,Vertex D,float * lambda,float * mu,int changed_AorB,int changed_C,int changed_D)4160 float _dxfSegmentIntersection(Vertex A, Vertex B, Vertex C, Vertex D,
4161 float *lambda, float *mu,
4162 int changed_AorB, int changed_C, int changed_D)
4163 {
4164 /* check whether for two segments in R3, the intersection is in the
4165 interior of each of them.
4166
4167 lamda values of 1 and 0 qualify as being in the interior, except
4168 when one of two edges is of length zero, in which case we set both
4169 lambda and mu to be -1, to reflect the fact that there is no
4170 intersection */
4171
4172 float dist = 0;
4173
4174 static VertexD w[3];
4175 static double work[2], max_z, min_z, sign = 0., scale;
4176
4177 register u_char i,j;
4178
4179
4180
4181 /* if both Vertices A and B did not change since the previous call
4182 of the routine, the scaling factor, householder vector and sign
4183 are still correct and do not need to be re-evaluated (this spares
4184 us the computation of one square root */
4185
4186 /* if C or D changed, we only have to recompute the portion
4187 corresponding to C or D
4188 in surface simplification, often the only variable that changes
4189 is D */
4190
4191 if (changed_AorB)
4192 {
4193
4194 for (j=0;j<3;j++)
4195 {
4196 /* substract A from the rest of the vertices */
4197
4198 w[0][j] = B[j]-A[j];
4199 w[1][j] = C[j]-A[j];
4200 w[2][j] = D[j]-A[j];
4201 }
4202
4203 /* record the norm of the first vector and scale the
4204 vertices accordingly */
4205
4206 scale = w[0][0] * w[0][0] + w[0][1] * w[0][1] + w[0][2] * w[0][2];
4207 if (scale > 0.0) scale = sqrt (scale);
4208
4209 /* particular case if the scale is zero */
4210
4211 if (scale > 0.)
4212 {
4213 /* apply a scaling factor so that w[0][2] will be negative */
4214
4215 sign = (w[0][2] > 0.) ? -1. : 1.;
4216
4217 sign /= scale;
4218
4219 for (i=0;i<3;i++)
4220 for (j=0;j<3;j++)
4221 w[i][j] *= sign; /* apply the scaling */
4222
4223 /* determine the Householder vector for transforming
4224 w[0] to z */
4225
4226 w[0][2] -= 1.;
4227
4228 _dxfHouseholderPreMultiplicationDbl((double *) w[1], 3,
4229 (double *) w[0], 3, 2, work);
4230
4231
4232 /* now the first segment is the same as (0,0,1) */
4233 }
4234 }
4235
4236 else
4237 {
4238 if (changed_C)
4239 {
4240 for (j=0;j<3;j++)
4241 w[1][j] = (C[j]-A[j]) * sign; /* apply scaling */
4242
4243 _dxfHouseholderPreMultiplicationDbl((double *) w[1], 3, /* apply Householder */
4244 (double *) w[0], 3, 1, work);
4245 }
4246
4247 if (changed_D)
4248 {
4249 for (j=0;j<3;j++)
4250 w[2][j] = (D[j]-A[j]) * sign; /* apply scaling */
4251
4252 _dxfHouseholderPreMultiplicationDbl((double *) w[2], 3, /* apply Householder */
4253 (double *) w[0], 3, 1, work);
4254 }
4255 }
4256
4257
4258 #define xC_ w[1][0]
4259 #define yC_ w[1][1]
4260 #define zC_ w[1][2]
4261 #define xD_ w[2][0]
4262 #define yD_ w[2][1]
4263 #define zD_ w[2][2]
4264
4265 if (scale > 0.)
4266 {
4267 /* there are two cases in which there is no intersection
4268 if zC and zD >1 or if zC and zD < 0 */
4269
4270 if (zC_> zD_) {max_z = zC_; min_z = zD_;}
4271 else {max_z = zD_; min_z = zC_;}
4272
4273 if ((min_z > 1.) || (max_z < 0.))
4274
4275 {
4276 *lambda = *mu = -1.;
4277 }
4278
4279 else
4280 {
4281 static double xDC, yDC, zDC, dlambda, dmu, a, half_b, c, tmp;
4282 /* x = xD + mu xDC
4283 y = yD + mu yDC */
4284
4285 xDC = xC_ - xD_ ;
4286 yDC = yC_ - yD_ ;
4287 zDC = zC_ - zD_ ;
4288
4289 c = xD_ * xD_ + yD_ * yD_;
4290 a = xDC * xDC + yDC * yDC;
4291
4292 if (a > 0.)
4293 {
4294 half_b = xD_ * xDC + yD_ * yDC;
4295
4296 dmu = -half_b / a;
4297 *mu = dmu; /* mu = -b/(2a) */
4298
4299 tmp = c + half_b * dmu;
4300 if (tmp >0.0) tmp = sqrt (tmp);
4301 else tmp= 0.0;
4302
4303 dist = scale * tmp; /* c - b^2/4a */
4304
4305 /* lambda is given by the 1-z_value at that point */
4306
4307 *lambda = 1. - zD_ - dmu * zDC;
4308 }
4309
4310 /* otherwise, either the two segment are parallel or the
4311 second segment is of zero length */
4312
4313 else if (min_z < max_z)
4314 {
4315
4316 /* the second segment is parallel to the first *and*
4317 there is an intersection */
4318
4319 min_z = MAX(min_z,0.);
4320 max_z = MIN(max_z,1.);
4321
4322 dlambda = (max_z + min_z)/2.;
4323
4324 /* this is definitely between [0,1] and corresponds to
4325 a point inside CD to */
4326
4327 *mu = (dlambda - zD_) / zDC;
4328
4329 *lambda = 1. - dlambda;
4330 dist = scale * sqrt (c);
4331
4332 }
4333 else
4334 {
4335 /* the second segment is degenerate */
4336 *lambda = *mu = -1.;
4337 }
4338
4339 }
4340 }/* endif (scale > 0) */
4341
4342 #undef xC_
4343 #undef yC_
4344 #undef zC_
4345 #undef xD_
4346 #undef yD_
4347 #undef zD_
4348
4349
4350 else
4351 {
4352 /* the first edge is degenerate. There is no mutual intersection */
4353
4354 *lambda = *mu = -1.;
4355
4356 }
4357
4358 return dist;
4359 }
4360
4361 /*+--------------------------------------------------------------------------+
4362 | |
4363 | _dxfHouseholderPreMultiplication |
4364 | |
4365 +--------------------------------------------------------------------------+*/
4366
_dxfHouseholderPreMultiplication(float * A,int mA,float * v,int m,int n,float * w)4367 int _dxfHouseholderPreMultiplication(
4368 float *A, /* input matrix, given column
4369 by column */
4370 int mA, /* leading dimension of matrix a */
4371 float *v, /* householder vector */
4372 int m, /* dimension of v */
4373 int n, /* second dimension of A */
4374 float *w /* workspace, dimension n */
4375 )
4376 {
4377 int retval = 0;
4378
4379 register int i,j ;
4380
4381 float *w1 = w;
4382
4383 /* 1) compute beta = -2 vT v: */
4384
4385 float beta = 0;
4386
4387 for (i=0;i<m;i++) beta += v[i]*v[i];
4388
4389 beta = -2./beta;
4390
4391 /* 2) compute w = beta AT v */
4392
4393 if (w1 == NULL)
4394 {
4395 w1 = (float *) DXAllocateZero (n * sizeof (float));
4396 }
4397
4398 if (w1)
4399 {
4400 float *A1 = A;
4401
4402 for (i=0;i< n; i++, A1 += mA)
4403
4404 {
4405 for (w1[i]=0., j=0; j<m;j++)
4406
4407 w1[i] += A1[j]*v[j];
4408
4409 w1[i] *= beta;
4410 }
4411
4412 /* 3) compute A = A + w vT */
4413
4414 A1 = A;
4415
4416 for (i=0;i< n; i++, A1 += mA)
4417
4418 for (j=0; j<m;j++)
4419
4420 A1[j] += w1[i] * v[j];
4421
4422 retval = 1;
4423
4424 if (w1 && w != w1)
4425 DXFree((Pointer) w1);
4426
4427 }
4428
4429 return retval;
4430 }
4431
4432 /*+--------------------------------------------------------------------------+
4433 | |
4434 | _dxfHouseholderPreMultiplicationDbl |
4435 | |
4436 +--------------------------------------------------------------------------+*/
4437
_dxfHouseholderPreMultiplicationDbl(double * A,int mA,double * v,int m,int n,double * w)4438 int _dxfHouseholderPreMultiplicationDbl(
4439 double *A, /* input matrix, given column
4440 by column */
4441 int mA, /* leading dimension of matrix a */
4442 double *v, /* householder vector */
4443 int m, /* dimension of v */
4444 int n, /* second dimension of A */
4445 double *w /* workspace, dimension n */
4446 )
4447 {
4448 int retval = 0;
4449
4450 register int i,j ;
4451
4452 double *w1 = w;
4453
4454 /* 1) compute beta = -2 vT v: */
4455
4456 double beta = 0;
4457
4458 for (i=0;i<m;i++) beta += v[i]*v[i];
4459
4460 beta = -2./beta;
4461
4462 /* 2) compute w = beta AT v */
4463
4464 if (w1 == NULL)
4465 {
4466 w1 = (double *) DXAllocateZero (n * sizeof (double));
4467 }
4468
4469 if (w1)
4470 {
4471 double *A1 = A;
4472
4473 for (i=0;i< n; i++, A1 += mA)
4474
4475 {
4476 for (w1[i]=0., j=0; j<m;j++)
4477
4478 w1[i] += A1[j]*v[j];
4479
4480 w1[i] *= beta;
4481 }
4482
4483 /* 3) compute A = A + w vT */
4484
4485 A1 = A;
4486
4487 for (i=0;i< n; i++, A1 += mA)
4488
4489 for (j=0; j<m;j++)
4490
4491 A1[j] += w1[i] * v[j];
4492
4493 retval = 1;
4494
4495 if (w1 && w != w1)
4496 DXFree((Pointer) w1);
4497
4498 }
4499
4500 return retval;
4501 }
4502
4503
4504 /*+--------------------------------------------------------------------------+
4505 | |
4506 | _dxfErrorWithinToleranceVBoundary2ndKnd |
4507 | |
4508 +--------------------------------------------------------------------------+*/
4509
_dxfErrorWithinToleranceVBoundary2ndKnd(SimpData * simp_data,int v0,int v1,int * vstar0,int * vstar1,int val0,int val1)4510 int _dxfErrorWithinToleranceVBoundary2ndKnd(SimpData *simp_data, int v0, int v1,
4511 int *vstar0, int *vstar1, int val0, int val1)
4512 {
4513 int edge_collapsible = 1;
4514
4515 /* we construct a tiling of (e*) and a tiling of (vs*)
4516 such that there is a one to one mapping between the two tilings */
4517
4518
4519 /* e*:
4520 v*1\ /v*0 \ /
4521 e*1\ *1/e*0 \ *1 /
4522 *val0-1 \ / *0 *0 \ / *val1-1
4523 ___________\/___________\/__________
4524 v*val0-1 v0 vs v1
4525 vl vr
4526 */
4527
4528 /* The difference between the interior vertex and boundary vertex
4529 approaches is that vs v0 and v1 are projected on edges exclusively.
4530 That limits the number of edge intersections that must be computed
4531 */
4532
4533 /* 1) determine the closest point to the vertex vs on the edge v0,v1 */
4534
4535 /* 2) determine the closest point to the vertex v0 on the edge
4536 (vl,vs) */
4537
4538 /* 3) determine the closest point to the vertex v1 on the edge (vs,vr) */
4539
4540 /* 4) take in turn each edge (vs,v*1) (vs,v*val0-1)
4541 and intersect them with each edge (v0,v*0), (v0, v*val0-1) */
4542
4543
4544 /* 5) do the same as 4) for the star of the vertex v1
4545 */
4546
4547 register int i;
4548
4549
4550 float
4551 potential_err_pivot,
4552 largest_err_pivot = 0.,
4553
4554 /* set the maximum tolerance for the star */
4555
4556 tol_pivot = MIN(simp_data->tol_volume[v0], simp_data->tol_volume[v1]);
4557
4558 for (i=0;i<val0;i++)
4559 if (simp_data->tol_volume[vstar0[i]] < tol_pivot)
4560 tol_pivot = simp_data->tol_volume[vstar0[i]];
4561
4562 for (i=1;i<val1;i++)
4563
4564 /* the first vertex of the star of v1 is the same as the first
4565 vertex of the star of v0: vstar1[0] = vstar0[0] */
4566
4567 if (simp_data->tol_volume[vstar1[i]] < tol_pivot)
4568 tol_pivot = simp_data->tol_volume[vstar1[i]];
4569
4570 /* verify that the maximum tolerance is not already exceeded at
4571 the vertices */
4572
4573 i = 0;
4574
4575 while (edge_collapsible && i < val0)
4576 /* the equality (star_err == tol) is allowed */
4577 if (simp_data->err_volume[vstar0[i++]] > tol_pivot)
4578 edge_collapsible = 0;
4579
4580 i = 1;
4581 while (edge_collapsible && i < val1)
4582 if (simp_data->err_volume[vstar1[i++]] > tol_pivot)
4583 edge_collapsible = 0;
4584
4585 if (simp_data->err_volume[v0] > tol_pivot || simp_data->err_volume[v1] > tol_pivot)
4586 edge_collapsible = 0;
4587
4588 if (edge_collapsible)
4589 {
4590 float
4591 t, dist;
4592
4593 Vertex wp[2];
4594
4595 bzero ((char *)wp, 2 *sizeof (Vertex));
4596
4597 _dxfErrorTestDataOnSurfaceInit(simp_data);
4598
4599 /* 1 determine the closest point to the vertex vs on the edge v0,v1 */
4600
4601 dist = _dxfClosestPointOnEdge(simp_data->simplified_vertex, wp[0],
4602 simp_data->vert[v0], simp_data->vert[v1],&t);
4603
4604 _dxfErrorTestDataOnSurfaceInterpolate2(simp_data, v0, v1, t);
4605
4606 potential_err_pivot = dist + t * simp_data->err_volume[v0] + (1.-t) * simp_data->err_volume[v1];
4607
4608 UPDATE_ERR_PIVOT;
4609
4610
4611 if (edge_collapsible)
4612 {
4613
4614 /* 2) determine the closest point to the vertex v0 on the edge
4615 (vl,vs) */
4616
4617
4618 dist = _dxfClosestPointOnEdge(simp_data->vert[v0], wp[0],
4619 simp_data->vert[vstar0[val0-1]], simp_data->simplified_vertex,
4620 &t);
4621
4622 /* t * err[vl] + 1-t * err[vs] >= dist + err[v0] */
4623
4624 if (t < 1.)
4625 {
4626 potential_err_pivot =
4627
4628 (dist + simp_data->err_volume[v0] -t * simp_data->err_volume[vstar0[val0-1]] )/ (1.-t);
4629
4630 UPDATE_ERR_PIVOT;
4631 }
4632
4633 /*
4634 v0
4635 if (1-t == 0 <==> t = 1) |
4636 that means that v0 projects orthogonally on vl vl____________vs
4637
4638 we can't treat this constraint and the edge is not collapsible
4639 unless err[vl] >= dist + err[v0]
4640 */
4641
4642 else if (simp_data->err_volume[vstar0[val0-1]] < dist + simp_data->err_volume[v0])
4643
4644 edge_collapsible = 0;
4645
4646 if (edge_collapsible)
4647 {
4648 if (simp_data->vx_data)
4649
4650 edge_collapsible =
4651 _dxfErrorTestDataOnSurfaceUpdate(simp_data,
4652 2, 1, vstar0[val0-1], 0, 1.-t, t, 0., /* vs, vl */
4653 v0, 0, 0, 1., 0., 0.);
4654
4655 /* 3) determine the closest point to the vertex v1 on the edge (vs,vr) */
4656
4657 dist = _dxfClosestPointOnEdge(simp_data->vert[v1], wp[0],
4658 simp_data->vert[vstar1[val1-1]], simp_data->simplified_vertex, &t);
4659
4660 /* t * err[vr] + 1-t * err[vs] >= dist + err[v1] */
4661
4662 if (t < 1.)
4663 {
4664 potential_err_pivot =
4665 (dist + simp_data->err_volume[v1] -t * simp_data->err_volume[vstar1[val1-1]])
4666 / (1.-t);
4667
4668 UPDATE_ERR_PIVOT;
4669 }
4670
4671 else if (simp_data->err_volume[vstar1[val1-1]] < dist + simp_data->err_volume[v1])
4672
4673 edge_collapsible = 0;
4674
4675
4676 if (simp_data->vx_data && edge_collapsible)
4677
4678 edge_collapsible = _dxfErrorTestDataOnSurfaceUpdate(simp_data,
4679 2, 1, vstar1[val1-1], 0, 1.-t, t, 0., /* vs, vl */
4680 v1, 0, 0, 1., 0., 0.);
4681
4682
4683 if (edge_collapsible)
4684 {
4685 float
4686 lambda1, lambda2;
4687
4688 int j;
4689
4690 /* 4) take in turn each edge (vs,v*1) (vs,v*val0-1)
4691 and intersect them with each edge (v0,v*0), (v0, v*val0-1) */
4692
4693 i = 1;
4694 while (edge_collapsible && i<val0)
4695 {
4696
4697 j = 0;
4698 while (edge_collapsible && j<val0)
4699 {
4700
4701 dist = _dxfSegmentIntersection(simp_data->simplified_vertex,
4702 simp_data->vert[vstar0[i]],
4703 simp_data->vert[v0],
4704 simp_data->vert[vstar0[j]],
4705 &lambda1, &lambda2, !j, 0, 1);
4706
4707 /* If the closest point is located in the
4708 inside of either edge, determine a crossing point. */
4709
4710 if (lambda2 >= 0. && lambda2 <= 1.) /* inside 2 */
4711 {
4712 dist += lambda2 * simp_data->err_volume[v0] +
4713 (1. -lambda2) * simp_data->err_volume[vstar0[j]];
4714
4715
4716 if (simp_data->vx_data && lambda1 >= 0. && lambda1 <= 1.)
4717
4718 edge_collapsible =
4719 _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 2,
4720 vstar0[i], 0, lambda1, 1. -lambda1, 0.,
4721 v0, vstar0[j], 0, lambda2, 1. -lambda2, 0.);
4722
4723
4724 if (lambda1 > 0. && lambda1 <= 1.) /* inside 1 */
4725
4726 {
4727 potential_err_pivot =
4728 (dist - (1. -lambda1) * simp_data->err_volume[vstar0[i]])
4729 / lambda1;
4730
4731
4732 UPDATE_ERR_PIVOT;
4733
4734 }
4735
4736 /* if lambda1 = 0 we can't accomodate
4737 that constraint with the error value at the simplified vertex
4738 and the edge might be declared "non-collapsible"
4739
4740 vs_______v*1
4741 |
4742 v0_|_v*0
4743
4744 dist == 0 can occur when the two vertices
4745 vstar0[j] and vstar0[i] are the same
4746 but we prohibit that */
4747
4748 else if (lambda1 == 0.)
4749 {
4750 /* the error at vs cannot accomodate
4751 that constraint, but maybe it is
4752 already satisfied */
4753
4754 if (simp_data->err_volume[vstar0[i]] < dist)
4755
4756 edge_collapsible = 0;
4757 }
4758 }
4759
4760
4761
4762 j++;
4763 if (i == j) j++; /* we prohibit (i == j) */
4764 }
4765 i++;
4766 }
4767
4768 if (edge_collapsible)
4769 {
4770 /* 5) do the same as 4) for the star of the vertex v1:
4771
4772 take in turn each edge (vs,v*1) (vs,v*val1-1)
4773 and intersect them with each edge (v1,v*0), (v1, v*val1-1) */
4774
4775 i = 1;
4776 while (edge_collapsible && i<val1)
4777 {
4778
4779 j = 0;
4780 while (edge_collapsible && j<val1)
4781 {
4782
4783 dist = _dxfSegmentIntersection(simp_data->simplified_vertex,
4784 simp_data->vert[vstar1[i]],
4785 simp_data->vert[v1],
4786 simp_data->vert[vstar1[j]],
4787 &lambda1, &lambda2, !j, 0, 1);
4788
4789
4790 if (lambda2 >= 0. && lambda2 <= 1.) /* inside 2 */
4791 {
4792
4793 dist +=
4794 lambda2 * simp_data->err_volume[v1] +
4795 (1. -lambda2) * simp_data->err_volume[vstar1[j]];
4796
4797
4798 if (simp_data->vx_data && lambda1 >= 0. && lambda1 <= 1.)
4799
4800 edge_collapsible =
4801 _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 2,
4802 vstar1[i], 0, lambda1, 1. -lambda1, 0.,
4803 v1, vstar1[j], 0, lambda2, 1. -lambda2, 0.);
4804
4805 if (lambda1 > 0. && lambda1 <= 1.) /* inside 1 */
4806
4807 {
4808 potential_err_pivot =
4809 (dist - (1. -lambda1) * simp_data->err_volume[vstar1[i]])
4810 / lambda1;
4811
4812
4813 UPDATE_ERR_PIVOT;
4814
4815 }
4816
4817 else if (lambda1 == 0.)
4818 {
4819 /* the error value at vs cannot accomodate that constraint,
4820 but maybe it is already satisfied */
4821
4822 if (simp_data->err_volume[vstar1[i]] < dist)
4823 edge_collapsible = 0;
4824 }
4825 }
4826
4827 j++;
4828 if (i == j) j++; /* we prohibit (i == j) */
4829 }
4830 i++;
4831 }
4832 }
4833 }
4834 }
4835 }
4836 }
4837
4838 /* update the tolerance and error volumes in case the edge will be simplified */
4839
4840 if (edge_collapsible)
4841 {
4842
4843 simp_data->err_volume[v0] = largest_err_pivot;
4844
4845 /* update the tolerance volume */
4846
4847 simp_data->tol_volume[v0] = MIN(simp_data->tol_volume[v0], simp_data->tol_volume[v1]);
4848
4849 _dxfErrorTestDataOnSurfaceEnd(simp_data, v0);
4850
4851 }
4852
4853 return edge_collapsible;
4854 }
4855
4856 /*+--------------------------------------------------------------------------+
4857 | |
4858 | _dxfRemoveEdgesFromHeap |
4859 | |
4860 +--------------------------------------------------------------------------+*/
4861
_dxfRemoveEdgesFromHeap(int * estar,u_short val,int * edge2index,SimpData * simp_data)4862 int _dxfRemoveEdgesFromHeap(int *estar, u_short val, int *edge2index,
4863 SimpData *simp_data)
4864 {
4865 int
4866 num_removed = 0;
4867
4868 register u_short i;
4869
4870 /* remove all edges of the star */
4871
4872 for(i=0;i<val;i++)
4873 {
4874 int the_edge = estar[i];
4875 int idx = edge2index[the_edge];
4876 num_removed += simp_data->remove_edge(simp_data, idx, the_edge);
4877 }
4878
4879 return num_removed;
4880 }
4881
4882 /*+--------------------------------------------------------------------------+
4883 | |
4884 | _dxfReinstateEdgesInHeap |
4885 | |
4886 +--------------------------------------------------------------------------+*/
4887
_dxfReinstateEdgesInHeap(int * estar,u_short val,SimpData * simp_data,int * num_edg_added)4888 int _dxfReinstateEdgesInHeap(int *estar, u_short val, SimpData *simp_data, int *num_edg_added)
4889
4890 /* return the number of edges that effectively have been reinstated
4891 (were in the queue before, were deleted and are now reintered */
4892 {
4893
4894 int
4895 num_edg_reinstated = 0;
4896
4897 register int i;
4898
4899 *num_edg_added = 0;
4900
4901 /* add all edges of the star to the heap */
4902
4903 for(i=0;i<val;i++)
4904 {
4905 int
4906 outside = (simp_data->edge2index[estar[i]] == simp_data->heap_outside_index),
4907 added = simp_data->add_edge(simp_data, estar[i]);
4908
4909 if (added)
4910 {
4911 *num_edg_added += 1;
4912 if (outside) num_edg_reinstated++;
4913
4914 }
4915 }
4916
4917 return num_edg_reinstated;
4918 }
4919
4920 /*+--------------------------------------------------------------------------+
4921 | |
4922 | _dxfCollapseBoundaryEdge2ndKnd |
4923 | |
4924 +--------------------------------------------------------------------------+*/
4925
_dxfCollapseBoundaryEdge2ndKnd(SimpData * simp_data,int edg_num,int v0,int v1,int val0,int val1,int * star0,int * star1,int * vstar0,int * vstar1,int * estar0,int * estar1,Vertex * s_normal0,Vertex * s_normal1,float * s_area0,float * s_area1)4926 int _dxfCollapseBoundaryEdge2ndKnd(SimpData *simp_data,
4927 int edg_num, int v0, int v1, int val0, int val1,
4928 int *star0, int *star1, int *vstar0, int *vstar1, int *estar0, int *estar1,
4929 Vertex *s_normal0, Vertex *s_normal1, float *s_area0, float *s_area1)
4930 {
4931 int
4932 n_edges_added[2],
4933 n_edges_reinstated = 0;
4934
4935
4936 /* e*:
4937 v*1\ /v*0 \ /
4938 e*1\ *1/e*0 \ *1 /
4939 *val0-1 \ / *0 *0 \ / *val1-1
4940 ___________\/___________\/__________
4941 v*val0-1 v0 vs v1
4942 vl vr
4943 */
4944
4945
4946 /* father relationships:
4947
4948 /\
4949 e*0[0] \e*1[0]
4950 / \
4951 /<-----\
4952 / ---\--> *1[1]
4953 *1[0]
4954 v0 <------- v1
4955
4956
4957 This is only true if val1 > 1
4958 otherwise if val1 ==1 we must have val0 > 1
4959 */
4960
4961 /* 1) update the vertex, triangle and edge fathers */
4962
4963 simp_data->vertex_father[v1] = v0;
4964
4965 if (val1 > 1)
4966 {
4967 simp_data->edge_father [estar1[0]] = simp_data->edge_father [estar0[0]];
4968
4969 simp_data->triangle_father[star1[0]] = simp_data->triangle_father[star1[1]];
4970 }
4971 else
4972 {
4973 simp_data->edge_father [estar0[0]] = simp_data->edge_father [estar1[0]];
4974
4975 simp_data->triangle_father[star0[0]] = simp_data->triangle_father[star0[1]];
4976 }
4977
4978 /* 2) change the valences of the vertices */
4979
4980 simp_data->valence[v0] = val0 + val1 - 2;
4981 simp_data->valence[vstar0[0]] --;
4982
4983 /* 3) remove from the edge heap all edges that have changed */
4984
4985 _dxfRemoveEdgesFromHeap(estar0, val0, simp_data->edge2index, simp_data);
4986
4987 _dxfRemoveEdgesFromHeap(estar1, val1, simp_data->edge2index, simp_data);
4988
4989 /* mark the 2 discarded edges.
4990
4991 verify that estar1[0] was not already counted as outside the heap */
4992
4993 if (simp_data->edge2index[estar1[0]] != simp_data->heap_outside_index)
4994 simp_data->num_edg_remaining_4_testing --;
4995
4996 simp_data->edge2index[estar1[0]] = -4;
4997 simp_data->edge2index[edg_num] = -7;
4998
4999
5000 /* 5) reinstate edges in the heap */
5001
5002 n_edges_reinstated =
5003 _dxfReinstateEdgesInHeap(estar0, val0, simp_data, n_edges_added);
5004
5005 /* reinstate all but the first edge (that collapses) in the star of v1 */
5006
5007 n_edges_reinstated +=
5008 _dxfReinstateEdgesInHeap(estar1+1, val1-1, simp_data, n_edges_added+1);
5009
5010 /* 6) update the number of edges left to be tested */
5011
5012 simp_data->num_edg_remaining_4_testing += n_edges_reinstated;
5013
5014 /* 7) update the coordinates of the simplified vertex */
5015
5016 memcpy(simp_data->vert[v0], simp_data->simplified_vertex, sizeof(Vertex));
5017
5018 /* replace the triangle normals, areas, and compactness in the simplified stars */
5019
5020 {
5021 register int i;
5022
5023 float
5024 *s_comp0 = s_area0 + val0,
5025 *s_comp1 = s_area1 + val1;
5026
5027 /* replace for v0 */
5028 for(i=1;i<val0;i++)
5029 {
5030 simp_data->area[star0[i]] = s_area0[i];
5031 simp_data->compactness[star0[i]] = s_comp0[i];
5032 memcpy(simp_data->normal[star0[i]], s_normal0[i], sizeof(Vertex));
5033 }
5034 /* replace for v1 */
5035
5036 for(i=1;i<val1;i++)
5037 {
5038 simp_data->area[star1[i]] = s_area1[i];
5039 simp_data->compactness[star1[i]] = s_comp1[i];
5040 memcpy(simp_data->normal[star1[i]], s_normal1[i], sizeof(Vertex));
5041 }
5042 }
5043
5044
5045
5046 return n_edges_added[0] + n_edges_added[1];
5047 }
5048
5049 /*+--------------------------------------------------------------------------+
5050 | |
5051 | _dxfCollapsibilityTestsBoundaryEdge2ndKind |
5052 | |
5053 +--------------------------------------------------------------------------+*/
5054
_dxfCollapsibilityTestsBoundaryEdge2ndKind(SimpData * simp_data,int edge_num,int v0,int v1,int val0,int val1,int move_simplified_vertex)5055 int _dxfCollapsibilityTestsBoundaryEdge2ndKind(SimpData *simp_data,
5056 int edge_num, int v0, int v1, int val0, int val1,
5057 int move_simplified_vertex)
5058 {
5059 int
5060 collapsible = 0,
5061 *edge_t = (int *) SurfEdgeGetTriangles (simp_data->edge_array+edge_num),
5062 *star0 = NULL;
5063
5064 Vertex *s_normal0 = NULL;
5065
5066 float *s_area0 = NULL;
5067
5068 /* perform collapsibility tests particular to a boundary edge of the
5069 second kind, such that both edge vertices are boundary vertices
5070
5071 In addition of the two vertices touching a boundary, the edge must
5072 effectively belong to the boundary, which means that exactly one
5073 triangle must be missing, with an index of -1. If both triangles had
5074 a -1 index, though, we would not have a valid edge */
5075
5076 if (simp_data->boundary_vert[v0] && simp_data->boundary_vert[v1] &&
5077 (edge_t[0] == -1 || edge_t[1] == -1))
5078 {
5079
5080 /* we only accept the edges such that both end vertices are located on the
5081 boundary */
5082
5083 int
5084 first_triangle = edge_t[0];
5085
5086 /* perform the valence test. For this particular case, the valence
5087 test is 1 <= v(v0) + v(v1) -2 <= valence_max */
5088
5089 if ((val0 + val1 -2 >= 1) && (val0 + val1 -2 < simp_data->valence_max))
5090 {
5091
5092 /* build the edge star:
5093 contrary to the procedure applied to the surface interior vertices,
5094 we compute the vertex star in a directed fashion.
5095 Here is the labeling system:
5096
5097 v*1\ /v*0 \ /
5098 e*1\ *1/e*0 \ *1 /
5099 *val0-1 \ / *0 *0 \ / *val1-1
5100 ___________\/___________\/__________
5101 v*val0-1 v0 v1
5102
5103 */
5104
5105 /* allocate storage space for the vertex stars */
5106
5107 star0 = (int *) simp_data->vertex_star_buffer;
5108
5109 bzero (star0, (val0 + val1) * 3 * sizeof (int));
5110
5111 /*
5112
5113 we next rotate the edge so that edge_t[0] is a valid triangle, t[0]
5114 i.e., so that the representation above holds true v[0]-----v[1]
5115 t[1]
5116 */
5117 if (first_triangle == -1)
5118 {
5119 /* rotate the edge, i.e. echange v0, v1, val0 and val1
5120 using first_triangle as a temporary storage location
5121
5122 assign edge_t[1] to first_triangle
5123 */
5124
5125 first_triangle = v0;
5126 v0 = v1;
5127 v1 = first_triangle;
5128
5129 first_triangle = val0;
5130 val0 = val1;
5131 val1 = first_triangle;
5132
5133 first_triangle = edge_t[1];
5134 }
5135
5136
5137 {int
5138 *vstar0 = star0 + val0,
5139 *estar0 = vstar0 + val0,
5140 *star1 = estar0 + val0,
5141 *vstar1 = star1 + val1,
5142 *estar1 = vstar1 + val1;
5143
5144 first_triangle = FATHER(first_triangle, simp_data->triangle_father);
5145
5146 _dxfDirectedParentVertexStar(v0, first_triangle, val0, star0, simp_data,
5147 DIRECTION_COUNTER_CLOCKWISE);
5148
5149 /* compute the star of v1 in reverse order */
5150
5151 _dxfDirectedParentVertexStar(v1, first_triangle, val1, star1, simp_data, DIRECTION_CLOCKWISE);
5152
5153
5154 /*
5155 up to now, we found a parent triangle
5156 and we located the next edge, "n_edgeCCW" or "n_edgeCW"
5157 depending on whether we are rotating counterclockwise
5158 or clockwise around "parent vertex".
5159
5160 The parent of this edge must be stored in the e_star.
5161
5162 _/
5163 n_edgeCCW/\
5164 parent vertex / \
5165 \_ /____\
5166 n_edge\_
5167 CW
5168
5169 Similarly to before, we need to orient the edge so that
5170 the parent of the first edge vertex is precisely "parent vertex"
5171 */
5172
5173 /* the manifold test is very similar to the test we
5174 perform for an interior edge case we compare the
5175 vertices of the links and determine whether duplicate
5176 vertices exist */
5177
5178 if (_dxfManifoldLinkTest(vstar0+1, vstar1+1, val0-1, val1-1))
5179 {
5180 Vertex
5181 *s_normal1;
5182 float
5183 *s_area1 ;
5184
5185 s_normal0 = (Vertex *) simp_data->edge_star_buffer_vx;
5186
5187 bzero (s_normal0, ( val0 + val1 ) * sizeof (Vertex));
5188
5189 s_normal1 = s_normal0 + val0;
5190
5191 s_area0 = (float *) simp_data->edge_star_buffer_fl;
5192
5193 bzero (s_area0, 2 * ( val0 + val1 ) * sizeof (float));
5194
5195 s_area1 = s_area0 + 2 * val0;
5196
5197 /* position the simplified vertex so that the length of the boundary
5198 will be preserved.
5199
5200 The vertex must be situated on an ellipse */
5201
5202 if ((_dxfSimplifiedBoundaryVertexLocation(simp_data, v0, v1, vstar0, vstar1,
5203 val0, val1, move_simplified_vertex) )
5204 &&
5205
5206 /* compute the difference in orientation for the
5207 triangle normals before and after simplification */
5208
5209 /* compute the minimum compactnesses before and after simplification */
5210
5211 (_dxfBoundaryCollapse2KndGeometricallyOk(simp_data, v0, v1, star0, vstar0, val0,
5212 s_normal0, s_area0,
5213 DIRECTION_COUNTER_CLOCKWISE,
5214 simp_data->min_scalprod_nor,
5215 simp_data->compactness_ratio) )
5216 &&
5217 (_dxfBoundaryCollapse2KndGeometricallyOk(simp_data, v1, v0, star1, vstar1, val1,
5218 s_normal1, s_area1, DIRECTION_CLOCKWISE,
5219 simp_data->min_scalprod_nor,
5220 simp_data->compactness_ratio) )
5221 ) /* geometry test */
5222
5223
5224 {
5225 /* compute the tolerance after simplification */
5226
5227 if (_dxfErrorWithinToleranceVBoundary2ndKnd(simp_data, v0, v1, vstar0, vstar1,
5228 val0, val1))
5229 {
5230
5231 simp_data->num_edg_collapsed +=2;
5232
5233 simp_data->num_edg_weights +=
5234
5235 _dxfCollapseBoundaryEdge2ndKnd(simp_data, edge_num, v0, v1, val0, val1,
5236 star0, star1, vstar0, vstar1,
5237 estar0, estar1, s_normal0, s_normal1,
5238 s_area0, s_area1);
5239
5240 collapsible = 1;
5241
5242 }
5243 else /* reject for tolerance off limits */
5244 simp_data->edg_rejected_4_tolerance ++;
5245
5246 }
5247 else /* reject for geometry */
5248 simp_data->edg_rejected_4_geometry ++;
5249 }
5250 else /* reject for topology */
5251
5252 {
5253 simp_data->edge2index[edge_num] = -3;
5254
5255 /* this condition will not be modified even if the
5256 surface is simplified around the edge */
5257
5258 simp_data->edg_rejected_4_topology ++;
5259 }
5260 }
5261 }
5262 else
5263 simp_data->edg_rejected_4_topology ++;
5264
5265 } /* end if (edge_t[0] == -1 || edge_t[1] == -1) */
5266
5267 return collapsible;
5268
5269 }
5270
5271 /*+--------------------------------------------------------------------------+
5272 | |
5273 | _dxfManifoldLinkTest1stKnd |
5274 | |
5275 +--------------------------------------------------------------------------+*/
5276
_dxfManifoldLinkTest1stKnd(int v0,int val0,int * star1,int * link1,int val1,SimpData * simp_data)5277 int _dxfManifoldLinkTest1stKnd(int v0, int val0, int *star1, int *link1, int val1, SimpData *simp_data)
5278 {
5279
5280 int link_is_manifold = 1;
5281
5282 /*
5283 v0 is a boundary vertex
5284 v1 is an interior vertex
5285
5286
5287 | /\vt
5288 |/ \<-----this triangle is *1[0]
5289 v0___\_v1
5290 |\ /
5291 | \ /<-----this triangle is *1[1]
5292 | vb
5293
5294
5295 this test works by taking
5296 all vertices in the link of v0, except v1, and verify that they are
5297 all different from the vertices of the link of v1, except vt, v0 and vb
5298
5299 we take from vstar0 as many vertices as val0 (remember that vstar0 has one element
5300 more than star0, because v0 is a boundary vertex)
5301
5302 we skip the vertices 0, 1, and val1-1 from the link of v1, i.e.,
5303 we take from vstar1 as many vertices as val1-3
5304
5305 */
5306
5307 int
5308 val = val0 + val1 -3;
5309
5310 /* since
5311 val0 + val1 <= SIMPLIFY_VALENCE_MAX4,
5312 I can use a static array here. */
5313
5314 static int
5315 link[SIMPLIFY_VALENCE_MAX4];
5316
5317 {
5318
5319 int
5320 i = 0,
5321
5322 tri_v_rotated[3], tri_v_parents[3],
5323
5324 /* rotated edge vertices and adjacent triangles: */
5325 edge_v_rotated[2], edge_v_parents[2], edge_t_rotated[2],
5326
5327 triangle = star1[0],
5328 edge;
5329
5330 do
5331 {
5332 triangle = FATHER (triangle, simp_data->triangle_father);
5333
5334 _dxfRotateParentTriangle(v0, simp_data->vertex_father,
5335 simp_data->triangles + 3 * triangle,
5336 tri_v_parents, tri_v_rotated);
5337
5338 edge = _dxfNextParentEdge(tri_v_rotated, DIRECTION_COUNTER_CLOCKWISE);
5339
5340 _dxfRotateParentEdge(v0, simp_data->vertex_father,
5341 (int *) SurfEdgeGetVertices (simp_data->edge_array + edge),
5342 (int *) SurfEdgeGetTriangles(simp_data->edge_array + edge),
5343 edge_v_parents, edge_v_rotated, edge_t_rotated);
5344
5345 link[i++] = edge_v_parents[1];
5346
5347 triangle = edge_t_rotated[0];
5348
5349 }
5350 while (triangle != -1);
5351
5352
5353 /* now repeat the same operation in the clockwise direction
5354 starting at vb (for Vertex_at_Bottom) */
5355
5356 triangle = star1[1];
5357
5358 do
5359 {
5360 triangle = FATHER (triangle, simp_data->triangle_father);
5361
5362 _dxfRotateParentTriangle(v0, simp_data->vertex_father,
5363 simp_data->triangles + 3 * triangle,
5364 tri_v_parents, tri_v_rotated);
5365
5366 edge = _dxfNextParentEdge(tri_v_rotated, DIRECTION_CLOCKWISE);
5367
5368 _dxfRotateParentEdge(v0, simp_data->vertex_father,
5369 (int *) SurfEdgeGetVertices (simp_data->edge_array + edge),
5370 (int *) SurfEdgeGetTriangles(simp_data->edge_array + edge),
5371 edge_v_parents, edge_v_rotated, edge_t_rotated);
5372
5373 link[i++] = edge_v_parents[1];
5374
5375 triangle = edge_t_rotated[1];
5376 }
5377 while (triangle != -1);
5378
5379
5380 /* gather vertices of link0 and link1 in a single array, sort
5381 the array, and determine whether there are any duplicate
5382 vertices */
5383
5384 memcpy(&link[i], link1+2, (val1 -3) * sizeof (int));
5385
5386 qsort(link, val, sizeof(int), (int (*)(const void*, const void*)) _dxfCmpIntSmallFirst);
5387
5388 i = 0;
5389
5390 while ((i < val -1 ) && (link_is_manifold))
5391 {
5392
5393
5394 if (link[i+1] == link[i]) /* if there are duplicate vertices,
5395 they must be contiguous in the
5396 sorted array */
5397
5398 link_is_manifold = 0;
5399 i++;
5400 }
5401
5402 }
5403
5404 return (link_is_manifold);
5405 }
5406
5407 /*+--------------------------------------------------------------------------+
5408 | |
5409 | _dxfBoundaryCollapse1stKndGeometricallyOk |
5410 | |
5411 +--------------------------------------------------------------------------+*/
5412
_dxfBoundaryCollapse1stKndGeometricallyOk(SimpData * simp_data,int v0,int * star0,int * vstar0,int val0,Vertex * s_normal0,float * s_area0,float min_scalprod_nor,float compactness_ratio)5413 int _dxfBoundaryCollapse1stKndGeometricallyOk(SimpData *simp_data, int v0, int *star0, int *vstar0,
5414 int val0, Vertex *s_normal0, float *s_area0,
5415 float min_scalprod_nor, float compactness_ratio)
5416
5417 {
5418 int collapsible = 1, i = 1;
5419
5420 Vertex faceVert[3];
5421
5422 float
5423 *s_comp0 = s_area0 + val0,
5424 min_compactness_before,
5425 min_compactness_after;
5426
5427 /*
5428 labeling system *[0]\ /
5429 v*[0]____\/____
5430 *[1]/\
5431 / \
5432 *[2]
5433
5434 both triangles *[0] and *[1] disappear after edge collapse
5435
5436 */
5437
5438
5439 /* compute the compactness of the collapsed triangles */
5440
5441 min_compactness_before = MIN ( simp_data->compactness[star0[1]], simp_data->compactness[star0[0]]);
5442
5443 min_compactness_after = 1; /* initialization of the compactness
5444 after simplification */
5445
5446 while (collapsible && i < val0-1)
5447 {
5448
5449 /* retrieve the face compactness before simplification */
5450 min_compactness_before = MIN ( simp_data->compactness[star0[i+1]], min_compactness_before);
5451
5452 /* recompute the compactness, and area, after changing the first vertex */
5453
5454 memcpy(faceVert[0], simp_data->simplified_vertex, sizeof (Vertex));
5455 memcpy(faceVert[1], simp_data->vert[vstar0[i]], sizeof (Vertex));
5456 memcpy(faceVert[2], simp_data->vert[vstar0[i+1]], sizeof (Vertex));
5457
5458 s_comp0[i+1] = _dxfFastCompactness3DTriangle2(faceVert, s_area0+i+1);
5459
5460 min_compactness_after = MIN ( s_comp0[i+1], min_compactness_after);
5461
5462 /* recompute the face normal assuming that the vertex v0 has
5463 moved and save this normal for future use*/
5464
5465 _dxfTriangleNormalQR2(faceVert, s_normal0[i+1]);
5466
5467 /* check whether the triangle normal orientation is consistent */
5468
5469 collapsible = (SCALPROD3(s_normal0[i+1], simp_data->normal[star0[i+1]]) > min_scalprod_nor);
5470
5471 i++;
5472 }
5473
5474 if (collapsible)
5475
5476 /* check whether the minimum triangle compactness is too much degraded */
5477
5478 collapsible = min_compactness_after > (compactness_ratio * min_compactness_before);
5479
5480 return collapsible;
5481 }
5482
5483 /*+--------------------------------------------------------------------------+
5484 | |
5485 | _dxfMakeBarycenter |
5486 | |
5487 +--------------------------------------------------------------------------+*/
5488
_dxfMakeBarycenter(int nV,Vertex * vv,Vertex bary,float * bary_coord)5489 void _dxfMakeBarycenter(int nV, Vertex *vv, Vertex bary, float *bary_coord)
5490 {
5491 register int i;
5492 register u_char j;
5493
5494 /* initialize barycenter */
5495
5496 for(j=0;j<3;j++)
5497 bary[j]=0.0;
5498
5499 for(i=0;i<nV;i++)
5500 {
5501 for(j=0;j<3;j++)
5502 bary[j]+= bary_coord[i]*vv[i][j];
5503 }
5504
5505 return;
5506 }
5507
5508 /*+--------------------------------------------------------------------------+
5509 | |
5510 | _dxfSolveSymmetric2x2Eqn |
5511 | |
5512 +--------------------------------------------------------------------------+*/
5513
_dxfSolveSymmetric2x2Eqn(double a,double b,double d,double e,double f,double * x,double * y)5514 int _dxfSolveSymmetric2x2Eqn(double a, double b, double d, double e, double f, double *x, double *y)
5515 {
5516 int retval = 0;
5517
5518 static double eps = 1e-13;
5519
5520 if (fabs(b) < eps)
5521 {
5522 *x = (fabs(a) < eps) ? 0. : e / a;
5523 *y = (fabs(d) < eps) ? 0. : f / d;
5524
5525 retval = 2;
5526 }
5527 else
5528 {
5529 double
5530 theta = (d - a) / 2. / b,
5531 sign = (theta > 0.) ? 1. : -1,
5532 t = sign / ( theta * sign + sqrt (theta * theta + 1.0)),
5533 /* we do not check for overflows of theta *theta */
5534 c = 1 / (sqrt (1.0 + t*t)),
5535 s = t * c,
5536
5537 /* | a b | | c s | | A | | c -s |
5538 | | = | | | | | |
5539 | b d | | -s c | | C | | s c |
5540 */
5541
5542 A = a - t * b,
5543 C = d + t * b,
5544 Rx = (fabs(A) < eps) ? 0. : (c * e - s * f) / A,
5545 Ry = (fabs(C) < eps) ? 0. : (s * e + c * f) / C;
5546
5547 *x = c * Rx + s * Ry ; /* apply the inverse rotation */
5548 *y = -s * Rx + c * Ry ; /* " " */
5549
5550 retval = 1;
5551 }
5552
5553 return retval;
5554 }
5555
5556 /*+--------------------------------------------------------------------------+
5557 | |
5558 | _dxfTriangle3DBarycentricCoordinates2 |
5559 | |
5560 +--------------------------------------------------------------------------+*/
5561
_dxfTriangle3DBarycentricCoordinates2(Vertex * tri,Vertex w,Vertex wp,Vertex bary,float * residual)5562 int _dxfTriangle3DBarycentricCoordinates2(Vertex *tri, Vertex w, Vertex wp, Vertex bary, float *residual)
5563
5564 /* compute the barycentric coordinates of a vector with respect to a point
5565 the projection wp and the residual are only valid if the function
5566 returns 1 (vertex projects inside)
5567 */
5568 {
5569 int inside;
5570
5571 static Vertex x1, x2;
5572
5573 static float eps_vertex_outside = EPS_VERTEX_OUTSIDE_TRIANGLE;
5574
5575 double a = 0., b = 0., c = 0., d = 0., e = 0., x3, x, y;
5576
5577 register u_char j;
5578
5579 int origin;
5580
5581 /*
5582 x1 is the vector representing the shortest edge of the triangle,
5583 which is best numerically
5584
5585 Min
5586
5587 || | x1[0] x2[0] | | w[0] - tri[origin][0] | ||^2
5588 || | x1[1] x2[1] | * | b1 | - | w[1] - tri[origin][1] | ||
5589 || | x1[2] x2[2] | | b2 | | w[2] - tri[origin][2] | ||
5590
5591
5592 using the normal equations
5593 */
5594
5595 _dxfTriangleBasisVectors(tri, x1, x2, origin);
5596
5597 for (j=0; j<3;j++)
5598 {
5599 x3 = w[j] - tri[origin][j];
5600
5601 a += x1[j] * x1[j];
5602 b += x1[j] * x2[j];
5603 c += x2[j] * x2[j];
5604 d += x1[j] * x3;
5605 e += x2[j] * x3;
5606 }
5607
5608 _dxfSolveSymmetric2x2Eqn(a,b,c,d,e, &x, &y);
5609
5610 bary [origin] = 1. - x - y;
5611
5612 bary [(origin + 1)%3] = x;
5613
5614 bary [(origin + 2)%3] = y;
5615
5616 _dxfMakeBarycenter(3, tri, wp, bary);
5617
5618 /* determine whether w projects inside or outside tri */
5619
5620 inside = !(bary[0] < eps_vertex_outside ||
5621 bary[1] < eps_vertex_outside ||
5622 bary[2] < eps_vertex_outside);
5623
5624 /* we actually do not need to create a variable eps_vertex_outside
5625 the following test would work all the same
5626
5627 inside = !(bary[0] < EPS_VERTEX_OUTSIDE_TRIANGLE ||
5628 bary[1] < EPS_VERTEX_OUTSIDE_TRIANGLE ||
5629 bary[2] < EPS_VERTEX_OUTSIDE_TRIANGLE);
5630
5631 */
5632
5633 *residual = DIST3(w, wp);
5634
5635 return inside;
5636
5637 }
5638
5639 /*+--------------------------------------------------------------------------+
5640 | |
5641 | _dxfClosestPointOnTriangle |
5642 | |
5643 +--------------------------------------------------------------------------+*/
5644
_dxfClosestPointOnTriangle(Vertex * tri,Vertex w,Vertex wp,Vertex bary,float * dist)5645 int _dxfClosestPointOnTriangle(Vertex *tri, Vertex w, Vertex wp, Vertex bary, float *dist)
5646 {
5647 int inside = _dxfTriangle3DBarycentricCoordinates2(tri, w, wp, bary, dist);
5648
5649 if (!inside)
5650 {
5651 static float eps= EPS_VERTEX_OUTSIDE_TRIANGLE;
5652
5653 /* either 2 barycentric coordinates are negative */
5654
5655 if (bary[0] < eps && bary[1] < eps)
5656 {
5657 bary[0] = bary[1] = 0.;
5658 bary[2] = 1.;
5659 memcpy(wp, tri[2], sizeof(Vertex));
5660 *dist = DIST3(w,wp);
5661 }
5662
5663 else if (bary[0] < eps && bary[2] < eps)
5664 {
5665 bary[0] = bary[2] = 0.;
5666 bary[1] = 1.;
5667 memcpy(wp, tri[1], sizeof(Vertex));
5668 *dist = DIST3(w,wp);
5669 }
5670
5671 else if (bary[1] < eps && bary[2] < eps)
5672 {
5673 bary[1] = bary[2] = 0.;
5674 bary[0] = 1.;
5675 memcpy(wp, tri[0], sizeof(Vertex));
5676 *dist = DIST3(w,wp);
5677 }
5678
5679 /* or 1 barycentric coordinate is negative */
5680
5681 else if (bary[0] < eps)
5682 {
5683 *dist = _dxfClosestPointOnEdge(w, wp, tri[1], tri[2], &bary[1]);
5684 bary[2] = 1. - bary[1];
5685 bary[0] = 0.;
5686 }
5687
5688 else if (bary[1] < eps)
5689 {
5690 *dist = _dxfClosestPointOnEdge(w, wp, tri[0], tri[2], &bary[0]);
5691 bary[2] = 1. - bary[0];
5692 bary[1] = 0.;
5693 }
5694
5695 else if (bary[2] < eps)
5696 {
5697 *dist = _dxfClosestPointOnEdge(w, wp, tri[0], tri[1], &bary[0]);
5698 bary[1] = 1. - bary[0];
5699 bary[2] = 0.;
5700 }
5701 }
5702
5703 return (inside);
5704 }
5705
5706 /*+--------------------------------------------------------------------------+
5707 | |
5708 | _dxfErrorWithinToleranceVBoundary1stKnd |
5709 | |
5710 +--------------------------------------------------------------------------+*/
5711
_dxfErrorWithinToleranceVBoundary1stKnd(SimpData * simp_data,int v0,int v1,int * vstar1,int val1)5712 int _dxfErrorWithinToleranceVBoundary1stKnd(SimpData *simp_data, int v0, int v1, int *vstar1, int val1)
5713 {
5714 int edge_collapsible = 1;
5715
5716 register int i;
5717
5718 float
5719 tol_pivot, potential_err_pivot, largest_err_pivot = 0.,
5720 dist,
5721 lambda1, lambda2;
5722
5723 /* set the maximum tolerance for the star */
5724
5725 tol_pivot = simp_data->tol_volume[v1];
5726
5727 for (i=0;i<val1;i++)
5728 if (simp_data->tol_volume[vstar1[i]] < tol_pivot)
5729 tol_pivot = simp_data->tol_volume[vstar1[i]];
5730
5731 /* verify that the maximum tolerance is not already exceeded at
5732 the vertices */
5733
5734 i = 0;
5735 while (edge_collapsible && i < val1)
5736 if (simp_data->err_volume[vstar1[i++]] > tol_pivot)
5737 edge_collapsible = 0;
5738
5739 if (simp_data->err_volume[v1] > tol_pivot)
5740 edge_collapsible = 0;
5741
5742 if (edge_collapsible)
5743 {
5744 Vertex
5745 triangle[3],
5746 vs[1],/* local copy of the coordinates of the simplified vertex */
5747 wp[2],
5748 alpha[SIMPLIFY_VALENCE_MAX1]; /* val0 + val1 -4 <= valence max <=> val0, val1 <= valence_max +1
5749 because val0 >=3 and val1 >=3 */
5750 float
5751 dp;
5752
5753 int
5754 i_min = -1;
5755
5756 bzero((char *)triangle, 3 * sizeof (Vertex));
5757
5758 bzero((char *)vs, sizeof (Vertex));
5759
5760 bzero((char *)wp, 2 * sizeof (Vertex));
5761
5762 bzero((char *)alpha, SIMPLIFY_VALENCE_MAX1 * sizeof (Vertex));
5763
5764 _dxfErrorTestDataOnSurfaceInit(simp_data);
5765
5766 /* the simplified is already positioned, since it must be v0 */
5767 memcpy(vs[0], simp_data->simplified_vertex, sizeof (Vertex));
5768
5769 /* if there is data attached to the surface, we also copy the data of v0
5770 to the new potential data */
5771
5772 if (simp_data->vx_data) memcpy(simp_data->vx_data_potential_values,
5773 simp_data->vx_data + v0*simp_data->data_dim,
5774 simp_data->data_dim * sizeof(float));
5775
5776 /* 1) determine the closest point to the vertex v1 on the triangles of star1 after simplification
5777 */
5778
5779 memcpy(triangle[0], vs[0], sizeof (Vertex));
5780
5781 /* v*[val1-1]
5782 labeling system *[0]\ /
5783 v*[0]____\/____
5784 *[1]/\
5785 v*[1]/ \
5786 *[2]
5787
5788 both triangles *[0] and *[1] disappear after edge collapse
5789
5790 */
5791 i = 1;
5792
5793 dist = 0.;
5794
5795 while (i<val1-1)
5796 {
5797 memcpy(triangle[1], simp_data->vert[vstar1[i]], sizeof(Vertex));
5798 memcpy(triangle[2], simp_data->vert[vstar1[i+1]], sizeof(Vertex));
5799
5800 _dxfClosestPointOnTriangle(triangle, simp_data->vert[v1], wp[0], alpha[i], &dp);
5801
5802 dp += simp_data->err_volume[v1];
5803
5804 if (alpha[i][0]>0.)
5805 {
5806 dp = (dp -
5807 alpha[i][1] * simp_data->err_volume[vstar1[i]] -
5808 alpha[i][2] * simp_data->err_volume[vstar1[i+1]])
5809 / alpha[i][0];
5810
5811 if (i_min == -1 || dp < dist) {dist = dp; i_min = i;}
5812
5813 }
5814 else
5815 {
5816 /* otherwise v1 projects onto the link edge of "triangle"
5817 the epsilon value at the simplified vertex cannot compensate
5818 but it is possible that the error values at the link are already sufficient */
5819
5820 if (alpha[i][1] * simp_data->err_volume[vstar1[i]] +
5821 alpha[i][2] * simp_data->err_volume[vstar1[i+1]] >= dp)
5822 {
5823 /* look no further for the closest triangle */
5824 i_min = i;
5825 i = val1-1; /* so that the loop will immediately stop*/
5826 dist = 0;
5827 }
5828 }
5829
5830 i++;
5831 }
5832
5833 /* the closest triangle wins */
5834
5835 if (i_min >=0)
5836
5837 {
5838 potential_err_pivot = dist;
5839
5840 UPDATE_ERR_PIVOT;
5841
5842 if (edge_collapsible && simp_data->vx_data)
5843
5844 edge_collapsible =
5845 _dxfErrorTestDataOnSurfaceUpdate(simp_data, 3, 1,
5846 vstar1[i_min], vstar1[i_min+1],
5847 alpha[i_min][0], alpha[i_min][1], alpha[i_min][2],
5848 v1, 0, 0, 1., 0., 0.);
5849
5850 }
5851
5852 else
5853 edge_collapsible = 0; /* v1 did not project inside of the new star */
5854
5855
5856 if (edge_collapsible)
5857 {
5858 /* 2)
5859
5860 / v*[val1-1] v*[val1-2]
5861 / __/ \ /
5862 / __/ \ /
5863 v0/__/________ ____\/____
5864 ==\ \__ v*0 /\
5865 vs \ \__ / \
5866 \ \ / \
5867 \ v*1 v*2
5868
5869 take in turn each edge (vs,v*2) ..., (vs, v*[val1-2])
5870
5871 and "intersect" them with each edge (v1, v*[1]) ... (v1, v*[val1-1])
5872 */
5873 int j;
5874
5875 i = 2;
5876
5877 while (edge_collapsible && i < val1-1)
5878 {
5879 j = 1;
5880
5881 while (edge_collapsible && j < val1)
5882 {
5883 dist = _dxfSegmentIntersection(simp_data->simplified_vertex,
5884 simp_data->vert[vstar1[i]],
5885 simp_data->vert[v1],
5886 simp_data->vert[vstar1[j]],
5887 &lambda1, &lambda2, (j==1), 0, 1);
5888
5889 /* If the closest point is located in the
5890 inside of either edge, determine a crossing point. */
5891
5892 if (lambda2 >= 0. && lambda2 <= 1.) /* inside 2 */
5893 {
5894 dist +=
5895 lambda2 * simp_data->err_volume[v1] +
5896 (1. -lambda2) * simp_data->err_volume[vstar1[j]];
5897
5898
5899 if (simp_data->vx_data && lambda1 >= 0. && lambda1 <= 1.)
5900
5901 edge_collapsible =
5902 _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 2,
5903 vstar1[i], 0, lambda1, 1. -lambda1, 0.,
5904 v1, vstar1[j], 0, lambda2, 1. -lambda2, 0.);
5905
5906
5907 if (lambda1 > 0. && lambda1 <= 1.) /* inside 1 */
5908
5909 {
5910 potential_err_pivot =
5911 (dist + (lambda1 -1.) * simp_data->err_volume[vstar1[i]])
5912 / lambda1;
5913
5914 UPDATE_ERR_PIVOT;
5915 }
5916 else if (lambda1 == 0.)
5917 {
5918 /* the error at vs cannot accomodate that constraint,
5919 but maybe it is already satisfied */
5920
5921 if (simp_data->err_volume[vstar1[i]] < dist)
5922 edge_collapsible = 0;
5923 }
5924 }
5925 j++;
5926
5927 if (i == j) j++; /* we prohibit (i == j) */
5928
5929 }
5930 i++;
5931 }
5932 }
5933
5934 }
5935
5936 if (edge_collapsible)
5937 {
5938 /* update the tolerance and error values */
5939
5940 /* only the error at the pivot will change */
5941
5942 simp_data->err_volume[v0] = largest_err_pivot;
5943
5944 /* update the tolerance volume */
5945
5946 simp_data->tol_volume[v0] = MIN(simp_data->tol_volume[v0], simp_data->tol_volume[v1]);
5947
5948 _dxfErrorTestDataOnSurfaceEnd(simp_data, v0);
5949
5950 }
5951
5952 return edge_collapsible;
5953 }
5954
5955 /*+--------------------------------------------------------------------------+
5956 | |
5957 | _dxfCollapseBoundaryEdge1stKnd |
5958 | |
5959 +--------------------------------------------------------------------------+*/
5960
_dxfCollapseBoundaryEdge1stKnd(SimpData * simp_data,int edg_num,int v0,int v1,int val0,int val1,int * star1,int * vstar1,int * estar1,Vertex * s_normal1,float * s_area1)5961 int _dxfCollapseBoundaryEdge1stKnd(SimpData *simp_data, int edg_num, int v0, int v1, int val0, int val1,
5962 int *star1, int *vstar1, int *estar1,
5963 Vertex *s_normal1, float *s_area1)
5964 {
5965 int
5966 n_edges_added = 0,
5967 n_edges_reinstated = 0,
5968 et, eb,
5969 tri_v_rotated[3], tri_v_parents[3];
5970
5971 /* v*[val1-1]
5972 labeling system *[0]\ /
5973 v*[0]____\/____
5974 *[1]/\
5975 v*[1]/ \
5976 *[2]
5977
5978 both triangles *[0] and *[1] disappear after edge collapse
5979
5980 */
5981
5982
5983 /* father relationships:
5984
5985 /\
5986 et / \e*1[val1-1]
5987 / \
5988 /<-----\
5989 / ---\--> *1[val1-1]
5990 *1[0]
5991 v0 <------- v1
5992 \ *1[1]/
5993 eb \ --/------> *1[2]
5994 \<---/e*1[1]
5995 \ /
5996 \/
5997 */
5998
5999 /* 1) update the vertex, triangle and edge fathers */
6000
6001 /* 1.1) First retrieve the edges et and eb on "top" and "bottom" of
6002 v0, since these edges were not pre-stored */
6003
6004 /* take the triangle star1[0], orient it with respect to v0
6005 and use the oriented triangle vertices to determine the next edge counter clockwise starting
6006 from v0: et */
6007
6008 _dxfRotateParentTriangle(v0, simp_data->vertex_father, simp_data->triangles + 3 * star1[0],
6009 tri_v_parents, tri_v_rotated);
6010
6011 et = (int) (((EdgeS *) _dxfFindEdge(simp_data->e_table, tri_v_rotated[0], tri_v_rotated[2])) -
6012 simp_data->edge_array);
6013
6014 et = FATHER (et, simp_data->edge_father);
6015
6016 /* do the same with the triangle star1[1] to determine the next edge clockwise starting
6017 from v0: eb */
6018
6019 _dxfRotateParentTriangle(v0, simp_data->vertex_father, simp_data->triangles + 3 * star1[1],
6020 tri_v_parents, tri_v_rotated);
6021
6022 eb = (int) (((EdgeS *) _dxfFindEdge(simp_data->e_table, tri_v_rotated[0], tri_v_rotated[1])) -
6023 simp_data->edge_array);
6024
6025 eb = FATHER (eb, simp_data->edge_father);
6026
6027 /* 1.2) */
6028
6029 simp_data->vertex_father[v1] = v0;
6030
6031 simp_data->edge_father [estar1[0]] = simp_data->edge_father [estar1[1]] = eb;
6032 simp_data->edge_father [estar1[val1-1]] = et;
6033
6034 simp_data->triangle_father[star1[0]] = star1[val1-1];
6035 simp_data->triangle_father[star1[1]] = star1[2];
6036
6037 /* 2) change the valences of the vertices */
6038
6039 simp_data->valence[v0] = val0 + val1 - 4;
6040 simp_data->valence[vstar1[1]] --;
6041 simp_data->valence[vstar1[val1-1]] --;
6042
6043 /* 3) remove from the edge heap all edges that have changed */
6044
6045 _dxfRemoveEdgesFromHeap(estar1, val1, simp_data->edge2index, simp_data);
6046
6047 /* mark the 2 discarded edges.
6048
6049 verify that they were not already counted as outside the heap */
6050
6051 if (simp_data->edge2index[estar1[1]] != simp_data->heap_outside_index)
6052 simp_data->num_edg_remaining_4_testing --;
6053
6054 if (simp_data->edge2index[estar1[val1-1]] != simp_data->heap_outside_index)
6055 simp_data->num_edg_remaining_4_testing --;
6056
6057 simp_data->edge2index[estar1[1]] = simp_data->edge2index[estar1[val1-1]] = -4;
6058 simp_data->edge2index[edg_num] = -7;
6059
6060
6061 /* 5) reinstate all but the two edges and last edge in the star of v1 */
6062
6063 n_edges_reinstated +=
6064 _dxfReinstateEdgesInHeap(estar1+2, val1-3, simp_data, &n_edges_added);
6065
6066 /* 6) update the number of edges left to be tested */
6067
6068 simp_data->num_edg_remaining_4_testing += n_edges_reinstated;
6069
6070 /* 7) update the coordinates of the simplified vertex" no need since
6071 its coordinates do not change*/
6072
6073 /* replace the triangle normals, areas and compactness in the simplified star of v1 */
6074
6075 {
6076 register int i;
6077
6078 float *s_comp1 = s_area1 + val1;
6079
6080 /* replace for v1 */
6081
6082 for(i=2;i<val1;i++)
6083 {
6084 simp_data->area[star1[i]] = s_area1[i];
6085 simp_data->compactness[star1[i]] = s_comp1[i];
6086 memcpy(simp_data->normal[star1[i]], s_normal1[i], sizeof(Vertex));
6087 }
6088 }
6089
6090 return n_edges_added;
6091 }
6092
6093
6094 /*+--------------------------------------------------------------------------+
6095 | |
6096 | _dxfCollapsibilityTestsBoundaryEdge1stKind |
6097 | |
6098 +--------------------------------------------------------------------------+*/
6099
_dxfCollapsibilityTestsBoundaryEdge1stKind(SimpData * simp_data,int edge_num,int v0,int v1,int val0,int val1)6100 int _dxfCollapsibilityTestsBoundaryEdge1stKind(SimpData *simp_data,
6101 int edge_num, int v0, int v1, int val0, int val1)
6102 {
6103 /* perform collapsibility tests particular to a boundary edge of
6104 the first kind, such that exactly one of the edge vertices is a
6105 boundary vertex */
6106
6107 int
6108 collapsible = 0;
6109
6110 int
6111 *star1 = NULL;
6112
6113 Vertex
6114 *s_normal1 = NULL;
6115
6116 float
6117 *s_area1 = NULL;
6118
6119 /* the valence of v0 becomes val0 + val1 - 4 */
6120
6121 if ((val0 + val1 -4 >= 1) && (val0 + val1 -4 < simp_data->valence_max))
6122 {
6123 int
6124 *edge_t = (int *) SurfEdgeGetTriangles (simp_data->edge_array+edge_num),
6125 first_tri = edge_t[0],
6126 second_tri = edge_t[1],
6127 /* this is the counter-clockwise order with respect to v1 */
6128 *vstar1,
6129 *estar1;
6130
6131 /* permute v0 and v1 so that v0 is the boundary vertex */
6132
6133 if (simp_data->boundary_vert[v1])
6134 {
6135 int
6136 tmp = v0;
6137 v0 = v1;
6138 v1 = tmp;
6139
6140 tmp = val0; val0 = val1; val1 = tmp;
6141 first_tri = edge_t[1]; second_tri = edge_t[0];
6142 }
6143
6144 /* initialize the storage pointer for the buffer star1 */
6145
6146 star1 = (int *) simp_data->vertex_star_buffer;
6147
6148 bzero (star1, 3 * val1 * sizeof (int));
6149
6150 /*
6151 | / first_triangle
6152 |/
6153 v0----------------v1
6154 \
6155 \
6156 boundary
6157 */
6158
6159
6160 estar1 = star1 + val1;
6161 vstar1 = estar1 + val1;
6162
6163 first_tri = FATHER(first_tri, simp_data->triangle_father);
6164
6165 /* first star elements (index 0)*/
6166 star1[0] = first_tri;
6167 estar1[0] = edge_num;
6168 vstar1[0] = v0;
6169
6170 /* complete the star from indices 1 to val1-1 */
6171 _dxfParentVertexStar(v1, second_tri, (u_short) val1, star1, simp_data);
6172
6173 if (_dxfManifoldLinkTest1stKnd(v0, val0, star1, vstar1, val1, simp_data))
6174 {
6175
6176 /* the simplified vertex is already positioned, since it must be v0 */
6177
6178 memcpy(simp_data->simplified_vertex, simp_data->vert[v0], sizeof (Vertex));
6179
6180 /* initialize the storage pointer for the buffer s_normal1 */
6181
6182 s_normal1 = (Vertex *) simp_data->edge_star_buffer_vx;
6183
6184 bzero(s_normal1, val1 * sizeof (Vertex));
6185
6186
6187 /* initialize the storage pointer for the buffer s_area1 */
6188
6189 s_area1 = (float *) simp_data->edge_star_buffer_fl;
6190
6191 bzero (s_area1, 2 * val1 * sizeof (float));
6192
6193
6194 if (_dxfBoundaryCollapse1stKndGeometricallyOk(simp_data, v1, star1, vstar1, val1,
6195 s_normal1, s_area1, simp_data->min_scalprod_nor,
6196 simp_data->compactness_ratio))
6197
6198 {
6199 /* perform the error tolerance test */
6200
6201 if (_dxfErrorWithinToleranceVBoundary1stKnd(simp_data, v0, v1, vstar1, val1))
6202 {
6203 simp_data->num_edg_collapsed +=3;
6204
6205 simp_data->num_edg_weights +=
6206
6207 _dxfCollapseBoundaryEdge1stKnd(simp_data, edge_num, v0, v1, val0, val1,
6208 star1, vstar1, estar1, s_normal1, s_area1);
6209
6210 collapsible = 1;
6211
6212 }
6213 else /* reject for tolerance off limits */
6214 simp_data->edg_rejected_4_tolerance ++;
6215 }
6216
6217 else /* reject for geometry */
6218 simp_data->edg_rejected_4_geometry ++;
6219
6220 }
6221 else
6222 {
6223 /* this condition will not change if the surface is
6224 further simplified */
6225
6226 simp_data->edge2index[edge_num] = -3;
6227
6228 simp_data->edg_rejected_4_topology ++;
6229 }
6230 }
6231 else
6232 simp_data->edg_rejected_4_topology ++;
6233
6234 return collapsible;
6235 }
6236
6237 /*+--------------------------------------------------------------------------+
6238 | |
6239 | _dxfCollapseTopologicallyFeasible |
6240 | |
6241 +--------------------------------------------------------------------------+*/
6242
_dxfCollapseTopologicallyFeasible(SimpData * simp_data,int * vstar0,int * vstar1,u_short val0,u_short val1)6243 int _dxfCollapseTopologicallyFeasible(SimpData *simp_data, int *vstar0, int *vstar1,
6244 u_short val0, u_short val1)
6245 {
6246
6247 register u_short i, j;
6248 int feasible = 1;
6249
6250 /* count each vertex only once
6251 *1[val1-1] /\ *0[1] counted twice
6252 *1[0] / \ *0[0]
6253 \ /
6254 *1[1] \/ *0[val0-1] counted twice
6255 */
6256
6257 {
6258 /* verify that no two vertices of the link are the same */
6259
6260 val0 -= 2;
6261 val1 -= 2;
6262
6263 i = 2;
6264 while(( i <= val0) && (feasible))
6265 {
6266
6267 j = 2;
6268 while ((j <= val1) && (feasible))
6269 {
6270
6271 if (vstar0[i] == vstar1[j])
6272 feasible = 0;
6273 j++;
6274 }
6275 i++;
6276 }
6277 }
6278
6279 return(feasible);
6280
6281 }
6282
6283 /*+--------------------------------------------------------------------------+
6284 | |
6285 | _dxfBuildParentEdgeStars |
6286 | |
6287 +--------------------------------------------------------------------------+*/
6288
_dxfBuildParentEdgeStars(EdgeS * edg,int v0,int v1,u_short val0,u_short val1,int * star0,int * star1,SimpData * simp_data)6289 void _dxfBuildParentEdgeStars(EdgeS *edg, int v0, int v1, u_short val0,
6290 u_short val1, int *star0, int *star1, SimpData *simp_data)
6291
6292
6293 /*...........................................................................*/
6294 /*
6295 The following procedure computes the star of v0 and v1 belonging to
6296 an "active" surface data structure, i.e. the surface is undergoing a
6297 simplification process.
6298
6299 the ordering of the stars is done such that:
6300
6301 e2 e1 e5 e4
6302 \ t2 / t1 \ t5 /
6303 t3 \ / t0 \ / t4
6304 e3 _______\/_______e0 e0 _______\/_______e3
6305 /\v0 v1 v0 /\v1
6306 t4 / \ t0 t1 / \ t3
6307 / t5 \ / t2 \
6308 e4 e5 e1 e2
6309 */
6310 /*...........................................................................*/
6311
6312 {
6313 int *estar0 = star0 + val0, *vstar0 = estar0 + val0;
6314 int *estar1 = star1 + val1, *vstar1 = estar1 + val1,
6315 e_num = (int) (edg - simp_data->edge_array);
6316
6317 int *edg_triangles = SurfEdgeGetTriangles(edg);
6318 /* complete first the straighforward portion of the stars
6319 /\
6320 / \
6321 / \ father(edg_triangles[0]);
6322 / \
6323 v0/edg->num\ v1
6324 \ /
6325 \ /
6326 \ / father(edg_triangles[1]);
6327 \ /
6328 \/
6329 */
6330
6331 star0 [0] = FATHER(edg_triangles[1], simp_data->triangle_father);
6332 vstar0[0] = v1;
6333 estar0[0] = e_num; /* the simplifiable edge is always a root edge */
6334
6335 _dxfParentVertexStar(v0, edg_triangles[0], val0, star0, simp_data);
6336
6337 star1 [0] = FATHER(edg_triangles[0], simp_data->triangle_father);
6338 vstar1[0] = v0;
6339 estar1[0] = e_num;
6340
6341 _dxfParentVertexStar(v1, edg_triangles[1], val1, star1, simp_data);
6342
6343 return;
6344 }
6345
6346 /*+--------------------------------------------------------------------------+
6347 | |
6348 | _dxfParentVertexStar |
6349 | |
6350 +--------------------------------------------------------------------------+*/
6351
_dxfParentVertexStar(int vf,int t0,u_short val,int * star,SimpData * simp_data)6352 void _dxfParentVertexStar(int vf, int t0, u_short val, int *star, SimpData *simp_data)
6353
6354 /*...........................................................................*/
6355 /*
6356 an "active" vertex is a vertex belonging to the inside of a surface
6357 ( not a boundary vertex ) undergoing a simplification process
6358 */
6359 /*...........................................................................*/
6360 {
6361 int *estar = star + val;
6362 int *vstar = estar + val;
6363 register int i;
6364 int t = t0;
6365
6366 /* complete the star from 1 to the valence,
6367 the 0-elements have been already completed */
6368
6369 for(i=1;i<val;i++)
6370 {
6371 register u_char j=0;
6372 int vA,vAf,vB,vBf,ef, e_num, *endpoints, *facing_tris;
6373 EdgeS *edg, *edgf;
6374
6375 t = FATHER(t, simp_data->triangle_father);
6376
6377 /* complete the triangle part of the vertex star */
6378 star[i] = t;
6379 do
6380 {
6381 vA = simp_data->triangles[3 * t + j];
6382 vB = simp_data->triangles[3 * t + (j+2)%3];
6383 vAf = FATHER(vA, simp_data->vertex_father);
6384 }
6385 while((j++ < 3) && (vAf != vf));
6386
6387 /* now complete the edge portion of the vertex star */
6388
6389 edg = _dxfFindEdge(simp_data->e_table,vA,vB);
6390
6391 e_num = (int) (edg - simp_data->edge_array); /* differential pointer is
6392 also the edge number */
6393
6394 ef = FATHER (e_num, simp_data->edge_father);
6395
6396 estar[i] = ef;
6397
6398 /* find out the parents of the vertices of ef */
6399
6400 edgf = simp_data->edge_array+ef;
6401
6402 endpoints = SurfEdgeGetVertices(edgf);
6403 facing_tris = SurfEdgeGetTriangles(edgf);
6404
6405 vAf = FATHER(endpoints[0], simp_data->vertex_father);
6406 vBf = FATHER(endpoints[1], simp_data->vertex_father);
6407
6408 /* and the vertex portion of the vertex star */
6409
6410 if (vAf == vf)
6411 {
6412 vstar[i] = vBf;
6413 /* /\f0 ^ccw direction
6414 the next triangle is v0/__\v1 | v0 = vf, v0 > v1
6415 \ / |
6416 \/f1 */
6417
6418 t = facing_tris[0];
6419 }
6420 else
6421 {
6422 vstar[i] = vAf;
6423 t = facing_tris[1];
6424 }
6425
6426 }
6427
6428 return;
6429 }
6430
6431 static void CentroidEdgeLink(Vertex centroid, Vertex *vert, int val0, int val02, int val014);
6432
6433 /*+--------------------------------------------------------------------------+
6434 | |
6435 | CentroidEdgeLink |
6436 | |
6437 +--------------------------------------------------------------------------+*/
6438
CentroidEdgeLink(Vertex centroid,Vertex * vert,int val0,int val02,int val014)6439 static void CentroidEdgeLink(Vertex centroid, Vertex *vert, int val0, int val02, int val014)
6440 {
6441 register int i;
6442 register u_char j;
6443
6444 bzero((char *)centroid, sizeof(Vertex));
6445
6446 for(j=0;j<3;j++)
6447 {
6448 for(i=0; i<=val02; i++)
6449
6450 centroid[j] += vert[i][j];
6451
6452 for(i=val0; i<=val014; i++)
6453
6454 centroid[j] += vert[i][j];
6455
6456 centroid[j] /= val014;
6457 }
6458
6459 return;
6460 }
6461
6462 static void BuildEdgeStarVal1Is3(Vertex *, Face *, int, int, int, int, int, int *, Vertex *);
6463
6464 /*+--------------------------------------------------------------------------+
6465 | |
6466 | BuildEdgeStarVal1Is3 |
6467 | |
6468 +--------------------------------------------------------------------------+*/
6469
BuildEdgeStarVal1Is3(Vertex * vert,Face * face,int v0,int v1,int val0,int val01,int val02,int * vstar0,Vertex * surfvert)6470 static void BuildEdgeStarVal1Is3(Vertex *vert, Face *face,
6471 int v0, int v1, int val0, int val01,
6472 int val02, int *vstar0, Vertex *surfvert)
6473 {
6474 register int i,i1;
6475
6476 /* example of numbering for the case when val0 = 6, val1 = 3
6477
6478 v1 _v0
6479 \ _//|
6480 \ t1 _/ / | The interior vertices come last
6481 t2 \ _/t0/ | The indices of triangles are shifted by one
6482 v2______\/___/ t6| with respect to the vertex star
6483 v5/\_ \v6 |
6484 t3 / \t5\ |
6485 / t4 \_\ |
6486 / \\|
6487 v3 v4 */
6488
6489
6490 /* 1) copy the locations of the vertices */
6491
6492 /* 1.1) star of v0 */
6493 for (i=0, i1=1; i< val01; i++, i1++)
6494 memcpy(vert[i], surfvert[vstar0[i1]], sizeof(Vertex));
6495 memcpy(vert[val01], surfvert[v0], sizeof(Vertex));
6496
6497 /* 1.2) star of v1 */
6498 memcpy(vert[val0], surfvert[v1], sizeof(Vertex));
6499
6500 /* 2) create the faces or connections between vertices */
6501
6502 face[0][0] = val01; face[0][1] = val0; face[0][2] = 0;
6503
6504 for (i = 1, i1 = 0; i <= val02; i++, i1++)
6505 {
6506 face[i][0] = val01; face[i][1] = i1; face[i][2] = i;
6507 }
6508 face[val01][0] = val01;
6509 face[val01][1] = val02;
6510 face[val01][2] = val0;
6511
6512 face[val0][0] = val0;
6513 face[val0][1] = val02;
6514 face[val0][2] = 0;
6515
6516 return;
6517 }
6518
6519
6520 /*+--------------------------------------------------------------------------+
6521 | |
6522 | _dxfBuildEdgeStar |
6523 | |
6524 +--------------------------------------------------------------------------+*/
6525
_dxfBuildEdgeStar(Vertex * star_vert,Face * star_face,Plane * star_plane,float * star_areas,float * star_comp,int v0,int v1,u_short val0,u_short val1,int * star0,int * star1,SimpData * s)6526 void _dxfBuildEdgeStar(Vertex *star_vert, Face *star_face, Plane *star_plane, float *star_areas,
6527 float *star_comp, int v0, int v1, u_short val0, u_short val1,
6528 int *star0, int *star1, SimpData *s)
6529
6530 /*...........................................................................*/
6531 /*
6532
6533 represent the edge star with origin at the centroid of the edge link
6534
6535 copy the vertices and faces of both vertex stars in a single array:
6536 here is the new numbering of the vertices:
6537
6538 1 0 val0+val1-4
6539 \ / \ / ..
6540 \ / \ /
6541 .. ____\/____ ____\/____val0+1
6542 val0-1/\v0 val0+val1-3/\
6543 / \ / \
6544 / \ / \val0
6545 .. val0 -2
6546
6547 The number of triangles in star_val is val0 + val1 -2
6548 The number of vertices is val0 + val1 -4 + 2,
6549 since v0 and v1 belong to the edge star
6550 In a vertex star, there is one more vertex (the center) than
6551 there are triangles.
6552 */
6553 /*...........................................................................*/
6554
6555 {
6556 /* there are two cases: either val1 is exactly 3 or it is > 3 */
6557
6558 int star_val = val0 + val1 - 2,
6559 val013 = star_val -1,
6560 val014 = val013 -1,
6561 val01 = val0 -1,
6562 val02 = val01 -1,
6563 *vstar0 = star0 + 2 * val0,
6564 *vstar1 = star1 + 2 * val1;
6565
6566 Vertex translation;
6567
6568 if (val1 == 3)
6569 {
6570 BuildEdgeStarVal1Is3(star_vert, star_face, v0, v1, val0, val01, val02, vstar0, s->vert);
6571 }
6572 else
6573 {
6574
6575 register int i,i1;
6576
6577 /* 0) copy the locations of the vertices */
6578
6579 /* 0.1) star of v0 */
6580
6581 for (i=0, i1=1; i< val01; i++, i1++)
6582 memcpy(star_vert[i], s->vert[vstar0[i1]], sizeof(Vertex));
6583
6584 memcpy(star_vert[val01], s->vert[v0], sizeof(Vertex));
6585
6586 /* 0.2) star of v1 */
6587
6588 for (i=val0, i1=2; i<= val014; i++, i1++)
6589 memcpy(star_vert[i], s->vert[vstar1[i1]], sizeof(Vertex));
6590
6591 memcpy(star_vert[val013], s->vert[v1], sizeof(Vertex));
6592
6593 /* 1) copy the faces */
6594
6595 /* numbering of the faces :
6596 val0 +
6597 \ 1 / \val1/
6598 2 \ / 0 \-3/
6599 ____\/____ ____\/____
6600 /\ /\
6601 .. / \ val0-1 / \val0+1
6602 / \ /val0\ */
6603
6604 star_face[0][0] = val01; star_face[0][1] = val013; star_face[0][2] = 0;
6605
6606 for (i = 1, i1 = 0; i <= val02; i++, i1++)
6607 {
6608 star_face[i][0] = val01; star_face[i][1] = i1; star_face[i][2] = i;
6609 }
6610 star_face[val01][0] = val01;
6611 star_face[val01][1] = val02;
6612 star_face[val01][2] = val013;
6613
6614 star_face[val0][0] = val013;
6615 star_face[val0][1] = val02;
6616 star_face[val0][2] = val0;
6617
6618
6619 for (i = val0+1, i1 = val0; i <= val014; i++, i1++)
6620 {
6621 star_face[i][0] = val013; star_face[i][1] = i1; star_face[i][2] = i;
6622 }
6623
6624 star_face[val013][0] = val013;
6625 star_face[val013][1] = val014;
6626 star_face[val013][2] = 0;
6627
6628 }
6629 /* 2) copy the normals of the faces and plug them into the
6630 equations of the planes */
6631
6632 _dxfCopyFaceNormalsAreasCompactness(star_plane, star_areas, star_comp,
6633 star0, star1, val0, val013, s->normal, s->area, s->compactness);
6634
6635 /* 3) compute the centroid of the edge link, and put it in the
6636 end of the star_vert array (the vertices v0 and v1 are excluded) */
6637
6638 CentroidEdgeLink(star_vert[star_val], star_vert, val0, val02, val014);
6639
6640 /* 4) change the coordinate system such that the origin will be the
6641 centroid of the edge link */
6642
6643
6644 _dxfOppositeVector(star_vert[star_val], translation);
6645
6646 _dxfApplyTranslation(star_val, star_vert, translation);
6647
6648 return;
6649 }
6650
6651 /*+--------------------------------------------------------------------------+
6652 | |
6653 | _dxfCopyFaceNormalsAreasCompactness |
6654 | |
6655 +--------------------------------------------------------------------------+*/
6656
_dxfCopyFaceNormalsAreasCompactness(Plane * plane,float * s_area,float * s_comp,int * star0,int * star1,int val0,int val013,Vertex * t_normal,float * t_area,float * t_comp)6657 void _dxfCopyFaceNormalsAreasCompactness(Plane *plane, float *s_area, float *s_comp,
6658 int *star0, int *star1, int val0, int val013,
6659 Vertex *t_normal, float *t_area, float *t_comp)
6660 {
6661 register int i, i1;
6662 int val01 = val0-1;
6663
6664 /* copy the normals of the faces and plug them into the
6665 equations of the planes */
6666
6667 for (i=0, i1=1; i< val01; i++, i1++)
6668 {
6669 memcpy(plane[i], t_normal[star0[i1]], sizeof(Vertex));
6670 s_area[i] = t_area[star0[i1]];
6671 s_comp[i] = t_comp[star0[i1]];
6672 }
6673
6674 memcpy(plane[val01], t_normal[star0[0]], sizeof(Vertex));
6675 s_area[val01] = t_area[star0[0]];
6676 s_comp[val01] = t_comp[star0[0]];
6677
6678 for (i=val0, i1=2; i<= val013; i++, i1++)
6679 {
6680 memcpy(plane[i], t_normal[star1[i1]], sizeof(Vertex));
6681 s_area[i] = t_area[star1[i1]];
6682 s_comp[i] = t_comp[star1[i1]];
6683 }
6684
6685 return;
6686 }
6687
6688 /*+--------------------------------------------------------------------------+
6689 | |
6690 | _dxfApplyTranslation |
6691 | |
6692 +--------------------------------------------------------------------------+*/
6693
_dxfApplyTranslation(int nV,Vertex * vv,Vertex T)6694 void _dxfApplyTranslation(int nV, Vertex *vv, Vertex T)
6695 {
6696 register int i;
6697 register u_char j;
6698
6699 for(i=0;i<nV;i++)
6700
6701 for(j=0;j<3;j++)
6702
6703 vv[i][j] = vv[i][j] + T[j];
6704
6705 return;
6706 }
6707
6708 /*+--------------------------------------------------------------------------+
6709 | |
6710 | _dxfOppositeVector |
6711 | |
6712 +--------------------------------------------------------------------------+*/
6713
_dxfOppositeVector(Vertex u,Vertex v)6714 void _dxfOppositeVector(Vertex u, Vertex v)
6715 {
6716 register u_char j;
6717
6718 for(j=0;j<3;j++)
6719
6720 v[j] = -u[j];
6721
6722 return;
6723 }
6724
6725 /*+--------------------------------------------------------------------------+
6726 | |
6727 | _dxfMinFloatArray |
6728 | |
6729 +--------------------------------------------------------------------------+*/
6730
_dxfMinFloatArray(float * array,int n)6731 float _dxfMinFloatArray(float *array, int n)
6732 {
6733 float min = array[0];
6734
6735 register int i;
6736
6737 for(i=1;i<n;i++)
6738 if (array[i] < min)
6739 min = array[i];
6740
6741 return min;
6742 }
6743
6744 /*+--------------------------------------------------------------------------+
6745 | |
6746 | _dxfStarPlaneEquationsAndVolume |
6747 | |
6748 +--------------------------------------------------------------------------+*/
6749
_dxfStarPlaneEquationsAndVolume(Vertex * star_vert,Face * star_face,Plane * star_plane,float * star_areas,int star_val)6750 double _dxfStarPlaneEquationsAndVolume(Vertex *star_vert, Face *star_face, Plane *star_plane,
6751 float *star_areas, int star_val)
6752 /*...........................................................................*/
6753 /*
6754 compute the equations of the planes of the star, supposing that the
6755 normals are known and are stored as the 3 first float values of the
6756 equations of the planes.
6757 compute the volume of the star
6758 */
6759 /*...........................................................................*/
6760 {
6761 double star_volume_x_6 = 0.;
6762
6763 register u_short j;
6764
6765 for(j=0;j<star_val;j++)
6766 {
6767 /* determine the equation of the plane when the normal is known.
6768
6769 (n,v) + d = 0 is the equation, n known vector, d unknown value.
6770 (n,v1) + d1 = 0;
6771 (n,v2) + d2 = 0; Three data points to determine the optimum d.
6772 (n,v3) + d3 = 0;
6773 d = -[ (n,v1) + (n,v2) + (n,v3) ]/3.;
6774 */
6775
6776 register u_char k;
6777
6778 star_plane[j][3] = 0.;
6779
6780 for(k=0;k<3;k++)
6781
6782 star_plane[j][3] -= SCALPROD3(star_plane[j], star_vert[star_face[j][k]]);
6783
6784 star_plane[j][3] /= 3.;
6785
6786
6787 /* for volume computation, retrieve the face area from the global array
6788 the volume is equal to one third of the triangle area times the
6789 height of the origin with respect to the triangle */
6790
6791 star_volume_x_6 -= star_areas[j] * star_plane[j][3] / 3.;
6792 }
6793
6794 return (star_volume_x_6 * 6.);
6795 }
6796
6797 /*+--------------------------------------------------------------------------+
6798 | |
6799 | _dxfComposeVectors |
6800 | |
6801 +--------------------------------------------------------------------------+*/
6802
_dxfComposeVectors(Vertex u,Vertex v,float lambda,Vertex w)6803 void _dxfComposeVectors(Vertex u, Vertex v, float lambda, Vertex w)
6804
6805 /*...........................................................................*/
6806 /*
6807 w = u + lambda v
6808 */
6809 /*...........................................................................*/
6810
6811 {
6812 register u_char j;
6813
6814 for(j=0;j<3;j++)
6815
6816 w[j] = u[j] + lambda * v[j];
6817
6818 return;
6819 }
6820
6821 /*+--------------------------------------------------------------------------+
6822 | |
6823 | _dxfPositionVertexLSConstraints2 |
6824 | |
6825 +--------------------------------------------------------------------------+*/
6826
_dxfPositionVertexLSConstraints2(Plane * star_plane,u_short star_val,Vertex simplified_vertex)6827 int _dxfPositionVertexLSConstraints2(Plane *star_plane, u_short star_val, Vertex simplified_vertex)
6828
6829 /*...........................................................................*/
6830 /*
6831
6832 Position the vertex such that
6833 1) the volume of the star is preserved
6834 2) the sum of distances plane-vertex is minimal plus
6835 3) the distance vertex-origin (minimize the sum of square distances
6836 vertex - vertices_of_the_link)
6837
6838 This version uses the Normal Least Squares Equations which have a
6839 closed form solution
6840
6841 We have changed the coordinate system such that the normal
6842 of the plane that guarantees to preserve the volume
6843 is given either by (0,0,1) or by (0,0,-1)
6844
6845 */
6846 /*...........................................................................*/
6847 {
6848 int retval = 1;
6849
6850 double
6851 A1[SIMPLIFY_VALENCE_MAX4],
6852 A2[SIMPLIFY_VALENCE_MAX4],
6853 B [SIMPLIFY_VALENCE_MAX4];
6854
6855 register int i;
6856
6857 /*1) each of the planes provide a LS equation, plus the minimum
6858 distance to the origin which provides 2 equations */
6859
6860 double
6861
6862 /* height of the plane (+/-) z + delta = 0 */
6863
6864 delta =
6865 star_plane[star_val][3] *
6866 star_plane[star_val][2];
6867
6868 /* 2) write the quadratic constraints associated with the planes */
6869
6870 for (i = 0; i < star_val; i++)
6871 {
6872 A1[i] = star_plane[i][0];
6873 A2[i] = star_plane[i][1];
6874 B [i] = delta * star_plane[i][2] - star_plane[i][3] ;
6875 }
6876
6877 /* 3) write the quadratic constraints associated with
6878 the center of the link which is the current origin */
6879
6880 A1[star_val ] = 1; /* x = 0 is the quadratic constraint */
6881 A1[star_val + 1] = 0;
6882 A2[star_val ] = 0;
6883 A2[star_val + 1] = 1; /* y = 0 is the quadratic constraint */
6884 B[star_val ] =
6885 B[star_val+1] = 0.;
6886
6887 {
6888 /* write the normal equations */
6889 static double a,b,c,d,e,x,y;
6890
6891 for (a=b=c=d=e=0.,i=0;i<star_val+2 ;i++)
6892 {
6893 a += A1[i] * A1[i];
6894 b += A1[i] * A2[i];
6895 c += A2[i] * A2[i];
6896 d += A1[i] * B [i];
6897 e += A2[i] * B [i];
6898 }
6899
6900 /* solve the normal equations */
6901
6902 _dxfSolveSymmetric2x2Eqn(a,b,c,d,e, &x, &y);
6903
6904 simplified_vertex[0] = (float) x;
6905
6906 simplified_vertex[1] = (float) y;
6907
6908 simplified_vertex[2] = (float) -delta;
6909
6910 }
6911
6912 return retval;
6913 }
6914
6915 /*+--------------------------------------------------------------------------+
6916 | |
6917 | _dxfFPlaneForIdenticalVolume |
6918 | |
6919 +--------------------------------------------------------------------------+*/
6920
_dxfFPlaneForIdenticalVolume(u_short val0,u_short val1,Plane p,Vertex * star_vert,Face * star_face,float star_volume)6921 void _dxfFPlaneForIdenticalVolume(u_short val0, u_short val1, Plane p, Vertex *star_vert,
6922 Face *star_face, float star_volume)
6923 {
6924 static Vertex n;
6925
6926 u_short i, val02 = val0 - 2, val013 = val0 + val1 - 3;
6927
6928 bzero((char *)p, sizeof(Vertex));
6929
6930 for(i=1; i<=val02; i++)
6931 {
6932 u_short v1 = star_face[i][1],
6933 v2 = star_face[i][2];
6934 VecProd3(star_vert[v1], star_vert[v2], n);
6935 AddVec(p, n);
6936 }
6937
6938 for(i=val0; i<=val013; i++)
6939 {
6940 u_short v1 = star_face[i][1],
6941 v2 = star_face[i][2];
6942 VecProd3(star_vert[v1], star_vert[v2], n);
6943 AddVec(p, n);
6944 }
6945
6946 p[3] = -star_volume;
6947
6948 /* normalize the equation of the plane */
6949 {
6950 double norm = NORM3(p);
6951
6952 if (norm > 0.0)
6953 for (i=0; i < 4; i++)
6954 p[i] /= norm;
6955 }
6956
6957 return;
6958 }
6959
6960
6961 /*+--------------------------------------------------------------------------+
6962 | |
6963 | _dxfSimplifiedVertexLocation |
6964 | |
6965 +--------------------------------------------------------------------------+*/
6966
_dxfSimplifiedVertexLocation(SimpData * s,u_short val0,u_short val1,u_short star_val,Vertex * star_vert,Face * star_face,Plane * star_plane,float * star_areas,float * work,u_char method)6967 void _dxfSimplifiedVertexLocation(
6968 SimpData *s,
6969 u_short val0,
6970 u_short val1,
6971 u_short star_val,
6972 Vertex *star_vert,
6973 Face *star_face,
6974 Plane *star_plane,
6975 float *star_areas,
6976
6977 /* workspace needed for the householder transform: */
6978
6979 float *work,
6980 u_char method)
6981 /*
6982 various methods for computing the position of the simplified vertex
6983 */
6984 {
6985 if (method == 0)
6986
6987 /* v0 is the simplified vertex */
6988 memcpy(s->simplified_vertex, star_vert[val0-1], sizeof(Vertex));
6989
6990
6991 else if (method == 1)
6992 {
6993 /* take the vertex that minimizes the sum of distances to the
6994 vertices of the edge star plus the planes defined by the
6995 triangles of the edge star and that preserves the volume of
6996 the edge star */
6997
6998 float
6999 star_volume, /* volume of the edge star (x 6) */
7000 sign;
7001
7002 int
7003 star_val_1 = star_val + 1,
7004 size_planes = (star_val_1) * sizeof (Plane);
7005
7006 float *star_plane_sav = (float *) (star_plane + star_val_1);
7007
7008 static Vertex
7009 v, z = {0,0,1.};
7010
7011 star_volume = _dxfStarPlaneEquationsAndVolume(star_vert, star_face, star_plane, star_areas, star_val);
7012
7013 /* compute the equation of the plane p such that the star of any
7014 vertex on that plane would have the same volume as star_volume */
7015
7016 _dxfFPlaneForIdenticalVolume(val0, val1, star_plane[star_val], star_vert, star_face, star_volume);
7017
7018 /* now change again the coordinate system such that either
7019 (0,0,1) or (0,0,-1) corresponds to the normal of P,
7020
7021 change the equations of the planes. THE VERTICES ARE LEFT UNCHANGED */
7022
7023 /* compute the householder vector */
7024
7025 sign = (star_plane[star_val][2] > 0.) ? 1. : -1.;
7026
7027 _dxfComposeVectors(star_plane[star_val], z, sign, v);
7028
7029 /* Apply the Householder transform to the planes
7030 while saving their previous positions */
7031
7032 memcpy(star_plane_sav, star_plane, size_planes);
7033
7034 _dxfHouseholderPreMultiplication((float *)star_plane, 4, v, 3, star_val_1, work);
7035
7036 /* this version of the positioning solves the Normal Equations
7037 rather than the LS equations and has a closed form solution
7038 (does not require usage of a numerical library) */
7039
7040 _dxfPositionVertexLSConstraints2( star_plane, star_val, s->simplified_vertex);
7041
7042 /* change back the position of the simplified vertex to
7043 the original coordinate system */
7044
7045 _dxfHouseholderPreMultiplication(s->simplified_vertex, 3, v, 3, 1, work);
7046
7047 /* idem with the planes */
7048
7049 memcpy(star_plane, star_plane_sav, size_planes);
7050
7051 }
7052
7053 else if (method == 2)
7054 {
7055 /* take the middle of v0 and v1 for the position of the
7056 simplified vertex */
7057 u_short val01 = val0-1, val013 = val0+val1-3;
7058
7059 register u_char j;
7060
7061 for (j=0;j<3;j++)
7062
7063 s->simplified_vertex[j] = (star_vert[val01 ][j] +
7064 star_vert[val013][j] ) /2.;
7065 }
7066
7067 return;
7068 }
7069
7070 /*+--------------------------------------------------------------------------+
7071 | |
7072 | _dxfCollapseGeometricallyFeasible |
7073 | |
7074 +--------------------------------------------------------------------------+*/
7075
_dxfCollapseGeometricallyFeasible(Vertex * svert,Face * sface,Plane * splane,float * old_comp,float new_comp_min,Vertex * snor0,Vertex * snor1,u_short val0,u_short sval,float min_scalprod,float compactness_ratio)7076 int _dxfCollapseGeometricallyFeasible(Vertex *svert, Face *sface,
7077 Plane *splane, float *old_comp, float new_comp_min,
7078 Vertex *snor0, Vertex *snor1,
7079 u_short val0, u_short sval,
7080 float min_scalprod, float compactness_ratio)
7081 {
7082 int feasible = 1;
7083
7084 /* check whether any triangle normal might have switched its orientation */
7085
7086 if (_dxfNoFlipOverCheck(1, val0-1, snor0, splane, min_scalprod) == 1 &&
7087 _dxfNoFlipOverCheck(val0, sval, snor1, splane, min_scalprod) == 1)
7088 {
7089 float old_comp_min = _dxfMinFloatArray(old_comp,sval);
7090
7091 feasible = (new_comp_min > (compactness_ratio * old_comp_min));
7092 }
7093 else feasible = 0;
7094
7095
7096 return feasible;
7097 }
7098
7099 /*+--------------------------------------------------------------------------+
7100 | |
7101 | _dxfNoFlipOverCheck |
7102 | |
7103 +--------------------------------------------------------------------------+*/
7104
_dxfNoFlipOverCheck(u_short first_face,u_short last_face,Vertex * nor,Plane * plane,float min_scalprod_nor)7105 int _dxfNoFlipOverCheck(u_short first_face, u_short last_face,
7106 Vertex *nor, Plane *plane, float min_scalprod_nor)
7107 {
7108
7109 register u_short i1=0, i=first_face;
7110
7111 int orientation_consistent;
7112
7113 do
7114 {
7115 orientation_consistent =
7116 (SCALPROD3(nor[i1],plane[i]) > min_scalprod_nor);
7117
7118 i++;
7119 i1++;
7120 }
7121 while ((orientation_consistent ) && (i < last_face));
7122
7123 return orientation_consistent;
7124
7125 }
7126
7127 /*+--------------------------------------------------------------------------+
7128 | |
7129 | _dxfErrorWithinToleranceV |
7130 | |
7131 +--------------------------------------------------------------------------+*/
7132
_dxfErrorWithinToleranceV(SimpData * simp_data,Vertex * star_vert,Face * star_face,Plane * star_plane,int * star0,int * star1,int val0,int val1)7133 int _dxfErrorWithinToleranceV(SimpData *simp_data, Vertex *star_vert,
7134 Face *star_face, Plane *star_plane,
7135 int *star0, int *star1, int val0, int val1)
7136
7137 {
7138 /* The purpose of Method V is to construct in parallel a tiling of
7139 e* (edge undergoing collapsibility tests) denoted til(e*) and a
7140 tiling of vs* (potential simplified vertex) denoted til(vs*), such
7141 that there is a one to one mapping between
7142 til(e*) and til(vs*)
7143 */
7144
7145 /* refresher: numbering strategies for e*:
7146
7147 numbering of the vertices :
7148 ~~~~~~~~~~~~~~~~~~~~~~~~~
7149
7150 1 0 val0+val1-4
7151 \ / \ / ..
7152 \ / \ /
7153 .. ____\/____ ____\/____val0+1
7154 val0-1 /\ val0+val1-3/\
7155 / \ / \
7156 / \ / \
7157 .. val0-2 val0
7158
7159 numbering of the faces : t0,t1,...t_val0+val1-3 triangles of e*
7160 ~~~~~~~~~~~~~~~~~~~~~~
7161 val0 +
7162 \ 1 / \val1/
7163 2 \ / 0 \-3/
7164 ____\/____ ____\/____
7165 /\ /\
7166 .. / \ val0-1 / \val0+1
7167 / \ /val0\
7168
7169 q1,...,qval0-2,....,qval0,...,q_val0+val1-3 triangles of vs*
7170 */
7171
7172 /* 1) determine the closest point to the vertex vs on the triangles t0
7173 and t_val0-1 */
7174
7175 /* 2) determine the closest point to the vertex v0 on the
7176 triangles q1,...,qval0-2. */
7177
7178 /* 3) determine the closest point to the vertex v1 on the triangles
7179 qval0,...,q_val0+val1-3 */
7180
7181 /* 4) Take in turn each edge (vs,v0),(vs,v1),..., (vs,v_val0-2)
7182 call the ith edge ei
7183 determine the closest point on each edge incident to v_val0-1:
7184
7185 (v_val0-1, v0),..., (v_val0-1, v_val0-2), and (v_val0-1, v_val0+val1-3)
7186 If the closest point is located in the inside of either edge,
7187 determine a crossing point.
7188 */
7189
7190
7191 /* 5) Take in turn each edge (vs,v_val0-2),...,(vs,v_val0 + val1-4) and
7192 (vs,v0),
7193
7194 call the ith edge ei
7195 determine the closest point on each edge incident to v_val0+val1-3:
7196
7197 */
7198
7199
7200 /* Firstly determine numbers of elements in old and new configurations */
7201
7202 int
7203 edge_collapsible = 1,
7204
7205 /* vertex stars */
7206 *vstar0 = star0 + 2 * val0,
7207 *vstar1 = star1 + 2 * val1,
7208 i,j, i_min = -1;
7209
7210 static Vertex
7211 vs[1], /* local copy of the simplified vertex */
7212 a_triangle[3],
7213 wp[2],
7214 alpha[SIMPLIFY_VALENCE_MAX1];
7215 /* val0 + val1 -4 <= valence max <=> val0, val1 <= valence_max +1
7216 because val0 >=3 and val1 >=3 */
7217
7218 static float
7219 dp[2];
7220
7221 float
7222 tol_pivot, potential_err_pivot, largest_err_pivot = 0., d_min,
7223 lambda1, lambda2;
7224
7225 _dxfErrorTestDataOnSurfaceInit(simp_data);
7226
7227 /* set the maximum tolerance for the star */
7228
7229 tol_pivot = simp_data->tol_volume[vstar0[0]];
7230
7231 for (i=1;i<val0;i++)
7232 if (simp_data->tol_volume[vstar0[i]] < tol_pivot)
7233 tol_pivot = simp_data->tol_volume[vstar0[i]];
7234
7235 for (i=0;i<val1;i++)
7236 if (simp_data->tol_volume[vstar1[i]] < tol_pivot)
7237 tol_pivot = simp_data->tol_volume[vstar1[i]];
7238
7239
7240 /* verify that the maximum tolerance is not already exceeded at
7241 the vertices */
7242
7243 i = 0;
7244 while (edge_collapsible && i < val0)
7245 /* the equality (star_err = tol) is allowed */
7246 if (simp_data->err_volume[vstar0[i++]] > tol_pivot)
7247 edge_collapsible = 0;
7248
7249 i = 0;
7250 while (edge_collapsible && i < val1)
7251 if (simp_data->err_volume[vstar1[i++]] > tol_pivot)
7252 edge_collapsible = 0;
7253
7254
7255 if (edge_collapsible)
7256 {
7257
7258 Vertex
7259 *qi, *t0, *t1;
7260
7261 /* local copy of the simplified vertex, translated by the due amount */
7262 memcpy(vs[0], simp_data->simplified_vertex, sizeof (Vertex));
7263 AddVec(vs[0], star_vert[val0+val1-2]);
7264
7265 /* 1) */
7266
7267 /* then fabricate the triangle t0 = v1 v0 vval0-1 */
7268 memcpy(a_triangle[0], simp_data->vert[vstar0[0]], sizeof (Vertex));
7269 t0 = ExtractTriangleFromVertexStar(vstar1,0);
7270
7271
7272 /* find the point closest to vs on t0 and record the projection as
7273 well as the distance */
7274
7275 _dxfClosestPointOnTriangle(t0, vs[0], wp[0], alpha[0], &dp[0]);
7276
7277 dp[0] +=
7278
7279 alpha[0][0] * simp_data->err_volume[vstar0[0]] +
7280 alpha[0][1] * simp_data->err_volume[vstar1[0]] +
7281 alpha[0][2] * simp_data->err_volume[vstar1[1]];
7282
7283 /* find the point closest to vs on t1 and record the same info
7284 triangle t1 = (v0, v1, v2) */
7285
7286 memcpy(a_triangle[0], simp_data->vert[vstar1[0]], sizeof (Vertex));
7287 t1 = ExtractTriangleFromVertexStar(vstar0,0);
7288
7289 _dxfClosestPointOnTriangle(t1, vs[0], wp[1], alpha[1], &dp[1]);
7290
7291 dp[1] +=
7292
7293 alpha[1][0] * simp_data->err_volume[vstar1[0]] +
7294 alpha[1][1] * simp_data->err_volume[vstar0[0]] +
7295 alpha[1][2] * simp_data->err_volume[vstar0[1]];
7296
7297 /* the triangle with the smallest potential error wins: */
7298
7299 if (dp[0] < dp[1])
7300
7301 {
7302 /* use the barycentric coordinates to interpolate the potential data values */
7303
7304 _dxfErrorTestDataOnSurfaceInterpolate3(simp_data,
7305 vstar0[0], vstar1[0], vstar1[1],
7306 alpha[0][0], alpha[0][1], alpha[0][2]);
7307
7308 potential_err_pivot = dp[0];
7309
7310 /* no need to test the data since we've just use this map point to
7311 compute new data values, so they have to match ! */
7312 /*
7313 if (simp_data->vx_data)
7314
7315 edge_collapsible = *//* SimpData simp_data, number_new_vertices_in_map, number_old_vertices,
7316 new_v2, new_v3, alpha_simp_vertex, alpha_n_v2, alpha_n_v3,
7317 old_v1, old_v2, old_v3, alpha_o_v1, alpha_o_v2, alpha_o_v3
7318 *//*
7319 _dxfErrorTestDataOnSurfaceUpdate(simp_data, 1, 3,
7320 0, 0, 1., 0., 0.,
7321 vstar0[0], vstar1[0], vstar1[1], alpha[0][0], alpha[0][1], alpha[0][2]);
7322 */
7323 }
7324 else
7325 {
7326 /* use the barycentric coordinates to interpolate the potential data values */
7327
7328 _dxfErrorTestDataOnSurfaceInterpolate3(simp_data,
7329 vstar1[0], vstar0[0], vstar0[1],
7330 alpha[1][0], alpha[1][1], alpha[1][2]);
7331
7332 potential_err_pivot = dp[1];
7333
7334 /* no need to test the data since we've just use this map point to
7335 compute new data values, so they have to match ! */
7336
7337 }
7338
7339 UPDATE_ERR_PIVOT;
7340
7341 /* 2) determine the closest point to the vertex v0 on the
7342 triangles q1,...,qval0-2. */
7343
7344 if (edge_collapsible)
7345 {
7346
7347
7348 /* the triangles q1,...,val0-2 are constructed explicitely
7349 with the triples of vertices:
7350 (vs,v0,v1), vs(v_val03, v_val02) */
7351
7352 memcpy(a_triangle[0], vs[0], sizeof(Vertex));
7353
7354 i = 1;
7355 d_min = 0;
7356
7357 while (i< val0-1)
7358 {
7359 qi = ExtractTriangleFromVertexStar(vstar0,i);
7360
7361 _dxfClosestPointOnTriangle(qi, simp_data->vert[vstar1[0]], wp[0], alpha[i], &dp[0]);
7362
7363 if (alpha[i][0]>0.)
7364
7365 /* If v0 does not project on the link, then the projection is useful and can produce
7366 a valid constraint */
7367
7368 {
7369 dp[0] =
7370 ( dp[0] + simp_data->err_volume[vstar1[0]]
7371 - alpha[i][1] * simp_data->err_volume[vstar0[i]]
7372 - alpha[i][2] * simp_data->err_volume[vstar0[i+1]])
7373 / alpha[i][0] ;
7374
7375 if (i_min == -1 || dp[0] < d_min) {d_min = dp[0]; i_min = i;}
7376 }
7377 else
7378 {
7379 /* otherwise v0 projects onto the link edge of "triangle"
7380 the epsilon value at the simplified vertex cannot compensate
7381 but it is possible that the error values at the link are high enough
7382 to treat that constraint */
7383
7384 if (alpha[i][1] * simp_data->err_volume[vstar0[i]] +
7385 alpha[i][2] * simp_data->err_volume[vstar0[i+1]] >= dp[0])
7386 {
7387 /* look no further for the closest triangle */
7388 i_min = i;
7389 i = val0-1; /* next i will be = val0 and the loop will immediately stop */
7390 d_min = 0.;
7391 }
7392 }
7393 i++;
7394 }
7395
7396 /* the triangle corresponding to the smallest error wins */
7397
7398 if (i_min >=0)
7399 {
7400 /* if vertex0 does not project on the rim,
7401
7402 i.e.,
7403
7404 if (alpha[i_min][0] > 0.)
7405
7406
7407 a potential new error value is generated at the pivot
7408 as follows: */
7409
7410
7411 potential_err_pivot = d_min;
7412
7413 UPDATE_ERR_PIVOT;
7414
7415 if (simp_data->vx_data && edge_collapsible)
7416
7417 edge_collapsible =
7418 _dxfErrorTestDataOnSurfaceUpdate(simp_data, 3, 1,
7419 vstar0[i_min], vstar0[i_min+1],
7420 alpha[i_min][0], alpha[i_min][1], alpha[i_min][2],
7421 vstar1[0], 0, 0, 1., 0., 0.);
7422
7423 /* even though v0 == vstar1[0] is the same vertex,
7424 the data is allowed to be modified at that vertex,
7425 we still have to perform this test in case the simplified
7426 vertex would have been allowed to move: we would then verify
7427 that the old data is still compatible with the new directions */
7428
7429 }
7430
7431 else edge_collapsible = 0; /* no projection was useful */
7432
7433 /* 4) */
7434
7435 if (edge_collapsible)
7436 {
7437 /* take each edge (vs,v1) .. (vs,v_val0-1)
7438 whereby v1, v_val0-1 refer to indices to
7439 the star the first vertex vstar0 */
7440
7441 i = 1;
7442 while (edge_collapsible && i<val0)
7443 {
7444
7445 j = 0;
7446 while (edge_collapsible && j<val0)
7447 {
7448
7449 d_min = _dxfSegmentIntersection(vs[0],
7450 simp_data->vert[vstar0[i]],
7451 simp_data->vert[vstar1[0]],
7452 simp_data->vert[vstar0[j]],
7453 &lambda1, &lambda2, !j, 0, 1);
7454
7455
7456 /* If the closest point is located in the
7457 inside of either edge, determine a crossing point. */
7458
7459 if (lambda2 >= 0. && lambda2 <= 1.) /* inside 2 */
7460 {
7461
7462
7463 d_min += lambda2 * simp_data->err_volume[vstar1[0]] +
7464 (1. -lambda2) * simp_data->err_volume[vstar0[j]];
7465
7466 if (simp_data->vx_data && lambda1 >= 0. && lambda1 <= 1.)
7467
7468 edge_collapsible =
7469
7470 _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 2,
7471 vstar0[i], 0, lambda1, 1. -lambda1, 0.,
7472 vstar1[0], vstar0[j], 0, lambda2, 1. -lambda2, 0.);
7473
7474 if (lambda1 > 0. && lambda1 <= 1.) /* inside 1 */
7475
7476 {
7477 potential_err_pivot =
7478 (d_min + ( lambda1 - 1.) * simp_data->err_volume[vstar0[i]])
7479 / lambda1;
7480
7481
7482 UPDATE_ERR_PIVOT;
7483
7484 }
7485
7486 /* LAMBDA1 = 0 IS AGAIN A SPECIAL CASE
7487 if lambda1 = 0
7488 we can't treat that constraint and the edge is
7489 not collapsible */
7490
7491 else if (lambda1 == 0.)
7492 {
7493 /* the error at vs cannot accomodate that constraint,
7494 but maybe it is already satisfied */
7495
7496 if (simp_data->err_volume[vstar0[i]] < d_min)
7497 edge_collapsible = 0;
7498 }
7499 }
7500
7501
7502
7503 j++; if (i==j) j++;
7504 }
7505 i++;
7506 }
7507
7508 /* 3) determine the closest point to the vertex v1 on
7509 the triangles qval0,...,q_val0+val1-3 */
7510
7511
7512 if (edge_collapsible)
7513 {
7514
7515 i_min = -1;
7516
7517 /* the triangles qval0,...,q_val0+val1-3 are
7518 constructed explicitely using the vertices of
7519 vstar1 and the triples of vertices with the
7520 triples of vertices: (vs,v1,v2), (vs,v_val1-2,
7521 v_val1-1) */
7522
7523 i = 1;
7524
7525 while (i<val1-1)
7526 {
7527
7528 qi = ExtractTriangleFromVertexStar(vstar1,i);
7529
7530 _dxfClosestPointOnTriangle(qi, simp_data->vert[vstar0[0]], wp[0], alpha[i], &dp[0]);
7531
7532 if (alpha[i][0]>0.)
7533
7534 {
7535 dp[0] =
7536 ( dp[0] + simp_data->err_volume[vstar0[0]]
7537 - alpha[i][1] * simp_data->err_volume[vstar1[i]]
7538 - alpha[i][2] * simp_data->err_volume[vstar1[i+1]])
7539 / alpha[i][0];
7540
7541 if (i_min == -1 || dp[0] < d_min) {d_min = dp[0]; i_min = i;}
7542 }
7543 else
7544 {
7545 /* otherwise v1 projects onto the link edge of "triangle"
7546 the epsilon value at the simplified vertex cannot compensate
7547 but it is possible that the error values at the link are high enough
7548 to treat that constraint */
7549 if (alpha[i][1] * simp_data->err_volume[vstar1[i]] +
7550 alpha[i][2] * simp_data->err_volume[vstar1[i+1]] >= dp[0])
7551 {
7552 /* look no further for the closest triangle */
7553 i_min = i;
7554 i = val1-1; /* and the loop will immediately stop */
7555 d_min = 0.;
7556 }
7557
7558 }
7559 i++;
7560 }
7561
7562 /* the closest triangle wins */
7563
7564 if (i_min >=0)
7565
7566 {
7567 potential_err_pivot = d_min;
7568
7569 UPDATE_ERR_PIVOT;
7570
7571
7572 if (simp_data->vx_data && edge_collapsible)
7573
7574 edge_collapsible =
7575 _dxfErrorTestDataOnSurfaceUpdate(simp_data, 3, 1,
7576 vstar1[i_min], vstar1[i_min+1],
7577 alpha[i_min][0], alpha[i_min][1], alpha[i_min][2],
7578 vstar0[0], 0, 0, 1., 0., 0.);
7579 }
7580 else edge_collapsible = 0; /* no projection was useful */
7581
7582 /* 5) */
7583
7584 if (edge_collapsible)
7585 {
7586 /* Take in turn each edge (vs,v2),(vs,v1),...,
7587 (vs,v_val1-2) wherein v2... refer to star 1
7588
7589 original vertices of the star:
7590 0,1,2,...,val1-1
7591 we take
7592 2,...,val1-2 */
7593
7594 i = 2;
7595
7596 while (edge_collapsible && i<val1-1)
7597 {
7598
7599 j = 0;
7600 while (edge_collapsible && j<val1)
7601 {
7602
7603 /* Determine the closest point on each
7604 edge incident to "v1" */
7605
7606 d_min = _dxfSegmentIntersection(vs[0],
7607 simp_data->vert[vstar1[i]],
7608 simp_data->vert[vstar0[0]],
7609 simp_data->vert[vstar1[j]],
7610 &lambda1, &lambda2, !j, 0, 1);
7611
7612
7613 if (lambda2 >= 0. && lambda2 <= 1.) /* inside 2 */
7614 {
7615 d_min += lambda2 * simp_data->err_volume[vstar0[0]] +
7616 (1.-lambda2) * simp_data->err_volume[vstar1[j]];
7617
7618 if (simp_data->vx_data && lambda1 >= 0. && lambda1 <= 1.)
7619
7620 edge_collapsible =
7621 _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 2,
7622 vstar1[i], 0, lambda1, 1. -lambda1, 0.,
7623 vstar0[0], vstar1[j], 0, lambda2, 1. -lambda2, 0.);
7624
7625 if (lambda1 > 0. && lambda1 <= 1.) /* inside 1 */
7626
7627 {
7628 potential_err_pivot =
7629 (d_min - (1.-lambda1) * simp_data->err_volume[vstar1[i]])
7630 / lambda1;
7631
7632
7633 UPDATE_ERR_PIVOT;
7634
7635 }
7636 else if (lambda1 == 0.)
7637 {
7638 /* the error at vs cannot accomodate that constraint,
7639 but maybe it is already satisfied */
7640
7641 if (simp_data->err_volume[vstar1[i]] < d_min)
7642 edge_collapsible = 0;
7643 }
7644 }
7645 j++;if (i==j) j++;
7646 }
7647 i++;
7648 }
7649 }
7650 }
7651 }
7652 }
7653 }
7654
7655
7656
7657 if (edge_collapsible)
7658 {
7659 /* update the tolerance and error values */
7660
7661 /* only the error at the pivot will change */
7662
7663 simp_data->err_volume[vstar1[0]] = largest_err_pivot;
7664
7665 /* update the tolerance volume */
7666
7667 simp_data->tol_volume[vstar1[0]] = MIN(simp_data->tol_volume[vstar1[0]],
7668 simp_data->tol_volume[vstar0[0]]);
7669
7670 _dxfErrorTestDataOnSurfaceEnd(simp_data, vstar1[0]);
7671
7672 }
7673
7674 return edge_collapsible;
7675 }
7676
7677 /*+--------------------------------------------------------------------------+
7678 | |
7679 | _dxfSimplifiedStarNormals |
7680 | |
7681 +--------------------------------------------------------------------------+*/
7682
_dxfSimplifiedStarNormals(Face * sface,Vertex * svert,Vertex * snor,float * sarea,float * scomp,Vertex simpvert,u_short first_face,u_short last_face)7683 void _dxfSimplifiedStarNormals(Face *sface, Vertex *svert, Vertex *snor,
7684 float *sarea, float *scomp, Vertex simpvert,
7685 u_short first_face, u_short last_face)
7686
7687 /*...........................................................................*/
7688 /*
7689 compute the face normals in the simplified star as well as the compactness
7690 and area of the faces.
7691 */
7692 /*...........................................................................*/
7693 {
7694 register u_short i,i1;
7695 register u_char j,k;
7696
7697 static VertexD triVert[3];
7698 static double n[4];
7699
7700 /* introduce the simplified vertex */
7701
7702 CopyType2Type(float, double, simpvert, triVert[0], 3);
7703
7704 for(i=first_face, i1=0; i < last_face; i++, i1++)
7705 {
7706 double sum_lengths_sq = 0.0;
7707
7708 /* get the second and third vertices of the face.. */
7709 u_short v1 = sface[i][1],
7710 v2 = sface[i][2];
7711
7712 /* ..and copy them at the 2nd and 3rd positions of faceVert */
7713
7714 CopyType2Type(float, double, svert[v1], triVert[1], 3);
7715
7716 CopyType2Type(float, double, svert[v2], triVert[2], 3);
7717
7718
7719 _dxfTriangleNormalQR2D(triVert, n);
7720
7721
7722 CopyType2Type(double, float, n, snor[i1], 3);
7723
7724 sarea[i1] = n[3];
7725
7726 /* compute the squared lenghts of the three edges: */
7727
7728 for(j=0;j<3;j++)
7729 {
7730 double len_edg_sq = 0, tmp;
7731
7732 for(k=0;k<3;k++)
7733 {
7734 tmp = triVert[(j+1)%3][k] - triVert[j][k];
7735 len_edg_sq += tmp * tmp;
7736 }
7737
7738 sum_lengths_sq += len_edg_sq;
7739
7740 }
7741
7742 scomp[i1] = FOUR_SQRT_3 * sarea[i1] / sum_lengths_sq;
7743
7744 }
7745
7746 return;
7747 }
7748
7749 /*+--------------------------------------------------------------------------+
7750 | |
7751 | _dxfReplaceFaceNormals |
7752 | |
7753 +--------------------------------------------------------------------------+*/
7754
_dxfReplaceFaceNormals(int * star,Vertex * new_nor,float * new_area,float * new_comp,Vertex * nor,float * area,float * comp,u_short val)7755 void _dxfReplaceFaceNormals(int *star, Vertex *new_nor, float *new_area, float *new_comp, Vertex *nor,
7756 float *area, float *comp, u_short val)
7757 {
7758
7759 /*...........................................................................*/
7760 /*
7761
7762 NOTE the count of the star goes from 0 to val-1
7763 the corresponding new normals go from 2 to val-1
7764 (triangle 0 and 1 are omitted since they will be removed)
7765 */
7766 /*...........................................................................*/
7767
7768 register u_short i,i1;
7769
7770 for(i=2, i1=0; i<val; i++, i1++)
7771 {
7772 memcpy(nor[star[i]], new_nor[i1], sizeof(Vertex));
7773 area[star[i]] = new_area[i1];
7774 comp[star[i]] = new_comp[i1];
7775 }
7776
7777 return;
7778 }
7779
7780 /*+--------------------------------------------------------------------------+
7781 | |
7782 | _dxfCollapseEdge |
7783 | |
7784 +--------------------------------------------------------------------------+*/
7785
_dxfCollapseEdge(SimpData * simp_data,int v0,int v1,int * star0,int * star1,u_short val0,u_short val1,Vertex * nor0,Vertex * nor1,float * area0,float * area1,float * comp0,float * comp1)7786 int _dxfCollapseEdge(SimpData *simp_data, int v0, int v1, int *star0, int *star1,
7787 u_short val0, u_short val1, Vertex *nor0, Vertex *nor1, float *area0, float *area1,
7788 float *comp0, float *comp1)
7789
7790 /*...........................................................................*/
7791 /*
7792
7793 Perform the simplification of an interior edge, and update the edge heap, by removing
7794 all edges of the star and adding back the new ones, keyed with the
7795 new distances.
7796
7797 father(*0[1]) /\
7798 /\/ \
7799 /\ \
7800 / *0[1]\
7801 /________\
7802 v0 \ /v1
7803 \ *O[0]/
7804 \ \/
7805 \ /\/
7806 \/ father(*0[1])
7807
7808
7809 / \
7810 The father of both / is \estar1[val1-1]
7811 /estar0[1] \
7812 / \
7813 /________ estar0[0] \
7814
7815
7816
7817 The father of / is \
7818 / \
7819 /estar1[1] \estar0[val0-1]
7820 / \
7821
7822 NOTE:
7823 In the *0 or star0 all triangles are parents of themselves
7824 (root of the disjoint-set forest)
7825
7826 in the surface, all vertices, triangles and edges are their original.
7827 We have to extract the parent each time.
7828 */
7829 /*...........................................................................*/
7830 {
7831 int
7832 *estar0 = star0 + val0, *vstar0 = estar0 + val0,
7833 *estar1 = star1 + val1, *vstar1 = estar1 + val1,
7834 n_edges_added[2],
7835 n_edges_reinstated = 0;
7836
7837 /* 1) perform the topological transformations */
7838
7839 simp_data->triangle_father[star1[0]] = star1[val1-1];
7840 simp_data->triangle_father[star1[1]] = star1[2];
7841
7842 simp_data->vertex_father[v1] = v0;
7843
7844 /* the remaining edges must be valid, and the way we choose to pick
7845 the triangle fathers dictates the way we choose edge father */
7846
7847 simp_data->edge_father[estar0[0]] =
7848 simp_data->edge_father[ estar1[val1-1] ] = estar0[1];
7849
7850 simp_data->edge_father[estar1[1]] = estar0[val0-1];
7851
7852 /* 2) change the valences of the vertices */
7853
7854 simp_data->valence[v0] = val0 + val1 - 4;
7855
7856 /* v0 is not the only vertex whose valence changes:
7857 the valence of vstar0[1] and vstar1[1] decreases by one
7858
7859 /\vstar0[1] |
7860 / \ |
7861 /____\ ---> |
7862 \ / |
7863 \ / |
7864 \/vstar1[1] | */
7865
7866 simp_data->valence[vstar0[1]] --;
7867 simp_data->valence[vstar1[1]] --;
7868
7869
7870 /* 3 remove the old edge lengths from the heap and put back the new ones */
7871
7872 /* 3.1 remove all edges in the star of v0 but v0 v1, because it was
7873 already removed from the heap (current edge under simplification) */
7874
7875 _dxfRemoveEdgesFromHeap(estar0+1, val0-1, simp_data->edge2index, simp_data);
7876
7877 /* 3.2 remove all edges but v0-v1 in the star of v1 */
7878
7879 _dxfRemoveEdgesFromHeap(estar1+1, val1-1, simp_data->edge2index, simp_data);
7880
7881 /* 4. mark the three edges that are to be discarded from the
7882 structure as collapsed with a -4 index */
7883
7884 /* update the number of edges left to be tested:
7885 the three collapsed edges wont be tested any more
7886 one among them was removed already : estar0[1]
7887 for the last two, we want to verify that they were not
7888 already counted as outside the heap */
7889
7890 if (simp_data->edge2index[estar1[val1-1]] != simp_data->heap_outside_index)
7891 simp_data->num_edg_remaining_4_testing --;
7892 if (simp_data->edge2index[estar1[1]] != simp_data->heap_outside_index)
7893 simp_data->num_edg_remaining_4_testing --;
7894
7895 simp_data->edge2index[estar1[val1-1]] = simp_data->edge2index[estar1[1]] = -4;
7896
7897 simp_data->edge2index[estar0[0]] = -7;
7898
7899 /* 5.3 reinstate in the heap all the edges but 0 in estar0
7900 that are their own parent, that is, all the edges that are a root
7901 of the disjoint set forest) */
7902
7903 n_edges_reinstated =
7904 _dxfReinstateEdgesInHeap(estar0+1, val0-1, simp_data, n_edges_added);
7905
7906 /* 2.3 reinstate in the heap all edges but 0 and 1 and val1-1 in estar1 */
7907
7908 n_edges_reinstated +=
7909 _dxfReinstateEdgesInHeap(estar1+2, val1-3, simp_data, n_edges_added+1);
7910
7911 /* update the number of edges left to be tested */
7912
7913 simp_data->num_edg_remaining_4_testing += n_edges_reinstated;
7914
7915 /* 3. copy the simplified vertex in the right place */
7916
7917 memcpy(simp_data->vert[v0], simp_data->simplified_vertex, sizeof(Vertex));
7918
7919 /* 4. replace the face normals, areas and compactness for the simplified stars with the new ones */
7920
7921 _dxfReplaceFaceNormals(star0, nor0, area0, comp0,
7922 simp_data->normal, simp_data->area, simp_data->compactness, val0);
7923 _dxfReplaceFaceNormals(star1, nor1, area1, comp1,
7924 simp_data->normal, simp_data->area, simp_data->compactness, val1);
7925
7926 return n_edges_added[0] + n_edges_added[1];
7927 }
7928
7929 /*+--------------------------------------------------------------------------+
7930 | |
7931 | _dxfSimplifyManifoldSurface |
7932 | |
7933 +--------------------------------------------------------------------------+*/
7934
_dxfSimplifyManifoldSurface(SimpData * simp_data)7935 int _dxfSimplifyManifoldSurface(SimpData *simp_data)
7936 {
7937 int
7938 e,
7939 num_edg_tested = 0, /* number of edges extracted from the heap */
7940 num_edg_added_to_heap = 0;
7941
7942 /* allocate buffers for storing various information on the vertex and edge stars */
7943
7944 /* since the valence of the simplified vertex can be as low as v0 + v1 - 4 (for an interior vertex)
7945 and is less than SIMPLIFY_VALENCE_MAX, assuming that v0 and v1 >= 1 we
7946 can be sure that the valence of v0 and v1 that are accepted for processing is <= SIMPLIFY_VALENCE_MAX+3.
7947 conservatively, we use SIMPLIFY_VALENCE_MAX+4 as a bound */
7948
7949 int
7950 vertex_star_buffer_size = 6 * SIMPLIFY_VALENCE_MAX4 * sizeof (int),
7951 edge_star_buffer_flsize = 9 * SIMPLIFY_VALENCE_MAX4 * sizeof (float),
7952 edge_star_buffer_vxsize = 4 * SIMPLIFY_VALENCE_MAX4 * sizeof (Vertex),
7953 edge_star_buffer_plsize = 4 * SIMPLIFY_VALENCE_MAX4 * sizeof (Plane),
7954 edge_star_buffer_fasize = 2 * SIMPLIFY_VALENCE_MAX4 * sizeof (Face),
7955
7956 *vertex_star_buffer = NULL;
7957
7958 Pointer
7959 edge_star_buffer_fl = NULL,
7960 edge_star_buffer_vx = NULL,
7961 edge_star_buffer_pl = NULL,
7962 edge_star_buffer_fa = NULL;
7963
7964 vertex_star_buffer = (int *) DXAllocateZero(vertex_star_buffer_size);
7965
7966 if (!vertex_star_buffer) goto error;
7967
7968 edge_star_buffer_fl = (Pointer) DXAllocateZero(edge_star_buffer_flsize);
7969
7970 if (!edge_star_buffer_fl) goto error;
7971
7972 edge_star_buffer_vx = (Pointer) DXAllocateZero(edge_star_buffer_vxsize);
7973
7974 if (!edge_star_buffer_vx) goto error;
7975
7976 edge_star_buffer_pl = (Pointer) DXAllocateZero(edge_star_buffer_plsize);
7977
7978 if (!edge_star_buffer_pl) goto error;
7979
7980 edge_star_buffer_fa = (Pointer) DXAllocateZero(edge_star_buffer_fasize);
7981
7982 if (!edge_star_buffer_fa) goto error;
7983
7984 /* reference such storage arrays in the simp_data structure */
7985
7986 simp_data->vertex_star_buffer = vertex_star_buffer;
7987 simp_data->edge_star_buffer_fl = edge_star_buffer_fl;
7988 simp_data->edge_star_buffer_vx = edge_star_buffer_vx;
7989
7990
7991 /* 1) mark the edges touching the surface boundary as not simplifiable: */
7992
7993 if (!simp_data -> simplify_boundary)
7994 /*num_edg_adjacent_2_boundary = */
7995 _dxfMarkEdgesAdjacent2Boundary(simp_data);
7996
7997 /* 2) fill a Heap with edges, keyed with the weight of the edge */
7998
7999
8000 while ((num_edg_added_to_heap=_dxfBuildEdgeHeap(simp_data))){
8001
8002 simp_data->num_edg_weights += num_edg_added_to_heap;
8003
8004 do
8005 {
8006
8007 e = simp_data->get_lowest_weight_edge(simp_data);
8008
8009 if (e != -1)
8010 {
8011
8012 /* variable set to 1 if the edge is effectively collapsed */
8013
8014 int
8015 collapsible = 0;
8016 EdgeS
8017 *edg = simp_data->edge_array+e;
8018 int
8019 *endpoints = SurfEdgeGetVertices(edg),
8020 v0 = FATHER(endpoints[0], simp_data->vertex_father),
8021 v1 = FATHER(endpoints[1], simp_data->vertex_father);
8022 u_short
8023 val0 = simp_data->valence[v0],
8024 val1 = simp_data->valence[v1];
8025
8026 num_edg_tested ++;
8027 /* mark the edge as removed from the heap: give an index such that HeapDelete would fail + 1*/
8028
8029 simp_data->edge2index[e] = simp_data->heap_outside_index;
8030
8031 /* decrement the number of edges remaining for testing */
8032
8033 simp_data->num_edg_remaining_4_testing --;
8034
8035 /* treat a boundary vertex of the second kind */
8036
8037 if (simp_data->simplify_boundary && simp_data->boundary_vert[v0] && simp_data->boundary_vert[v1])
8038
8039 collapsible =
8040 _dxfCollapsibilityTestsBoundaryEdge2ndKind(simp_data, e, v0, v1,
8041 val0, val1, simp_data->preserve_boundary_length);
8042
8043 /* treat a boundary vertex of the fisrt kind */
8044
8045 else if (simp_data->simplify_boundary && (simp_data->boundary_vert[v0] ||
8046 simp_data->boundary_vert[v1]))
8047
8048 collapsible =
8049 _dxfCollapsibilityTestsBoundaryEdge1stKind(simp_data, e, v0, v1, val0, val1);
8050
8051 /* treat an interior edge */
8052 else
8053 {
8054 /* perform the edge valence test for an interior edge */
8055
8056 if ( (val0 + val1 - 4 >= 3) && (val0 + val1 - 4 < simp_data->valence_max))
8057 {
8058
8059 /* BUILD THE 2 VERTEX STARS */
8060
8061 int
8062 *star0 = (int *) vertex_star_buffer,
8063 *estar0 = star0 + val0, *vstar0 = estar0 + val0,
8064 *star1 = vstar0 + val0, *estar1 = star1 + val1,
8065 *vstar1 = estar1 + val1;
8066
8067 bzero(star0, 3 * (val0 + val1) * sizeof (int));
8068
8069 /* 3) construct the edge star */
8070
8071 _dxfBuildParentEdgeStars(edg,v0,v1,val0,val1,star0,star1,simp_data);
8072
8073 /* 3.1) test whether the simplification is topologically feasible */
8074
8075 if (_dxfCollapseTopologicallyFeasible(simp_data, vstar0, vstar1, val0, val1))
8076 {
8077 /* 3.2) compute a simplified vertex, and test whether the simplification is
8078 geometrically possible, if not, go to END */
8079
8080 /* the number of vertices is given by val0 + val1 - 2 and the number of
8081 triangles is also val0 + val1 - 2 */
8082
8083 u_short
8084 sval = val0 + val1 -2,
8085 sval1 = sval + 1,
8086 val02 = val0-2, val12 = val1-2;
8087
8088 /* BUILD THE EDGE STAR */
8089
8090 Vertex
8091 *svert = (Vertex *) edge_star_buffer_vx,
8092 *snor0 = svert + sval1,
8093 *snor1 = snor0 + val02;
8094 Face
8095 *sface = (Face *) edge_star_buffer_fa;
8096 Plane
8097 *splane = (Plane *) edge_star_buffer_pl;
8098 float
8099 *scomp0 = (float *) edge_star_buffer_fl,
8100 *sarea0 = scomp0 + val02,
8101 *scomp1 = sarea0 + val02,
8102 *sarea1 = scomp1 + val12,
8103 *scomp = sarea1 + val12,
8104 *sarea = scomp + sval,
8105 *wspace = sarea + sval,
8106 min_comp_after;
8107
8108 bzero((char *)svert, edge_star_buffer_vxsize);
8109 bzero((char *)sface, edge_star_buffer_fasize);
8110 bzero((char *)splane, edge_star_buffer_plsize);
8111 bzero((char *)scomp0, edge_star_buffer_flsize);
8112
8113 /* the new origin is placed at the centroid of the edge link */
8114
8115 _dxfBuildEdgeStar(svert, sface, splane, sarea, scomp,
8116 v0, v1, val0, val1, star0, star1, simp_data);
8117
8118 _dxfSimplifiedVertexLocation(simp_data, val0, val1, sval, svert, sface,
8119 splane, sarea, wspace, simp_data->preserve_volume);
8120
8121 /* to determine whether the collapse is geometrically
8122 possible, the new vertex normals are needed */
8123
8124 _dxfSimplifiedStarNormals(sface, svert, snor0, sarea0, scomp0,
8125 simp_data->simplified_vertex, 1, val0-1);
8126 _dxfSimplifiedStarNormals(sface, svert, snor1, sarea1, scomp1,
8127 simp_data->simplified_vertex, val0, sval);
8128
8129 /* minimum triangle compactness after simplification :*/
8130 min_comp_after = MIN( _dxfMinFloatArray( scomp0, val02 ),
8131 _dxfMinFloatArray( scomp1, val12 ) );
8132
8133 if (_dxfCollapseGeometricallyFeasible(svert, sface, splane, scomp,
8134 min_comp_after, snor0, snor1, val0, sval,
8135 simp_data->min_scalprod_nor,
8136 simp_data->compactness_ratio))
8137
8138 {
8139 /* 3.3) if collapsible, test whether the new errors are within the
8140 tolerance limits */
8141
8142 collapsible = _dxfErrorWithinToleranceV(simp_data, svert, sface, splane,
8143 star0, star1, val0, val1);
8144
8145 if (collapsible)
8146
8147 {
8148
8149 /* 3.4) simplify (collapse the edge ) and replace old edges in
8150 the heap with new ones */
8151
8152 simp_data->num_edg_collapsed += 3;
8153
8154 /* translate back to the (0,0,0) origin */
8155
8156 AddVec( simp_data->simplified_vertex, svert[sval]);
8157
8158 simp_data->num_edg_weights +=
8159 _dxfCollapseEdge(simp_data, v0, v1, star0, star1,
8160 val0, val1, snor0, snor1, sarea0, sarea1,
8161 scomp0, scomp1);
8162 }
8163 else
8164 {
8165 /* rejected for tolerance off limits the edge can be
8166 reinserted afterwards if its length changes*/
8167
8168 simp_data->edg_rejected_4_tolerance ++;
8169 }
8170
8171 }
8172 else
8173 {/* geometrically infeasible */
8174
8175 simp_data->edg_rejected_4_geometry ++;
8176 }
8177
8178 }
8179 else
8180 { /* topologically infeasible */
8181
8182 simp_data->edge2index[e] = -3;
8183
8184 /* this condition will not be modified even if the surface is simplified around
8185 the edge */
8186
8187 simp_data->edg_rejected_4_topology ++;
8188 }
8189
8190
8191 }
8192 else
8193 {
8194
8195 /* If the valence is too high or too low, reject also for topology, However, this condition
8196 might change */
8197
8198 simp_data->edg_rejected_4_topology ++;
8199 }
8200 } /* end treat interior vertex */
8201
8202 } /* end if (e != -1) */
8203
8204 } while(e != -1);
8205
8206
8207 } /* end while (_dxfBuildEdgeHeap) */
8208
8209 DXFree((Pointer)vertex_star_buffer);
8210
8211 DXFree(edge_star_buffer_fl);
8212
8213 DXFree(edge_star_buffer_vx);
8214
8215 DXFree(edge_star_buffer_pl);
8216
8217 DXFree(edge_star_buffer_fa);
8218
8219
8220 return 1;
8221
8222 error:
8223
8224 if (vertex_star_buffer) DXFree((Pointer)vertex_star_buffer);
8225
8226 if (edge_star_buffer_fl) DXFree(edge_star_buffer_fl);
8227
8228 if (edge_star_buffer_vx) DXFree(edge_star_buffer_vx);
8229
8230 if (edge_star_buffer_pl) DXFree(edge_star_buffer_pl);
8231
8232 if (edge_star_buffer_fa) DXFree(edge_star_buffer_fa);
8233
8234 return 0;
8235 }
8236
8237 /*****************************************************************************/
8238
8239
8240