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 /*********************************************************************/ 13 /* 14 15 Author: Andre GUEZIEC, March/April 1997 16 17 */ 18 /*********************************************************************/ 19 20 21 #ifndef _SIMPLESURF_H_ 22 #define _SIMPLESURF_H_ 23 24 /* SimplifySurface 25 26 is a module that approximates a triangulated surface and resamples data 27 28 */ 29 30 /* partial bug log 31 ~~~~~~~~~~~~~~~ 32 33 this is a partial log of the bugs I discovered when developing this module. 34 it is there mainly to help me, and other potential developers working 35 on this module spend less time hunting bugs, especially if we want to 36 change the code and redo the same errors. 37 38 39 3/29/97 DXCull, DXInvalidateConnections, DXInvalidateUnreferencedPositions 40 41 I was trying to run something like: 42 43 obj1 = DXInvalidateUnreferencedPositions(obj) 44 obj2 = DXCull(obj1) 45 46 That corresponded to what the User manual suggested 47 but I was getting an error at DXGetObjectClass with the result of DXCull 48 49 I decided to code things differently: 50 51 new_in[0] = DXInvalidateConnections(in[0]); 52 if (!DXInvalidateUnreferencedPositions(new_in[0])) goto error; 53 if (!DXCull(new_in[0])) goto error; 54 55 56 ........................................................................... 57 58 59 3/31/97 DXNewArrayV must be used in place of DXNewArray if a "shape" int * vector 60 is passed to the routine rather than an int. Typically if rank > 1 but also if 61 rank = 1 and shape is obtained from DXArrayInfo, which expects an int *. 62 63 The problem with this bug is that the error message appeared in the end of the module 64 something like Out of Memory Out of Memory, and I had no obvious way of tracking it down to the 65 source. 66 67 68 ........................................................................... 69 70 71 3/31/97 _dxfVertexNormalsfromTriangleNormals(nor, new_nV, new_nT, new_connections, 72 face_normals, new_face_areas); 73 74 was executed *after* DXFree((Pointer)new_connections); which is a stupid thing to do! 75 I was freeing after 76 DXAddArrayData(con, 0, new_nT, new_connections); 77 and when later on I added the vertex normals from triangle normals routine, 78 I forgot that this array was freed. Part of the problem is that I wasnt used to the 79 style of allocating and freeing with DXAllocate and DXFree whereby you do a goto error 80 if something bad happens. I was doing things differently: 81 82 return code = 0 83 buffer = (malloc) 84 85 if (buffer) { do the work; free(buffer); set return code to 1} 86 else { write error message } 87 88 return(return code) 89 90 91 ........................................................................... 92 93 94 95 3/31/97 DXConvertArray() sets an error code if the array is already of the same type that requires 96 conversion. Same thing: it is not documented and it took me a while to realize what was going wrong 97 since the error was showing up in the end of the module execution, and was not informative: 98 99 Out of Memory, or Null Object, or something similar 100 101 102 ........................................................................... 103 104 105 106 4/2/97 MACRO DEFINITIONS: I had a hard time compiling (2 hours) because of syntax errors 107 in the macro definitions that are not picked up by the compiler. They are detected much 108 later and point so seemingly incomprehensible errors in the c-code. 109 110 111 ........................................................................... 112 113 114 4/2/97 SurfEdgeLabel: when building the "e_table" hash table of edges, I was forgetting 115 to initialize the label of the edge to 1 116 117 SurfEdgeLabel(edge) = 1; 118 119 Accordingly, edges that were boundary edges had a label of 28800 or something like that 120 (an arbitrary value) and were misclassified as interior edges. 121 122 then, the routine that was looping around edges was crashing because it expected boundary 123 edges to be interior edges. 124 125 This bug was easy to locate in dbx: I started from the location where the program crashed, 126 I printed the edge array and noticed right away that 127 the labels were funny. I was leaded quickly to the routine that generated the edges 128 129 130 131 ........................................................................... 132 133 134 4/2/97 DXAllocate the valence instead of DXAllocateZero: stupid error. I had only two 135 DXAllocate instead of DXAllocateZero in the entire program. 136 137 I noticed it in dbx because the module was simplifying surface so little. 138 by stopping at the right line, I notices way too many topologically infeasible collapses. 139 That made me want to check the valence 140 141 simp_data->boundary_vert was also DXAllocated and that was probably another bug, fixed 142 by hitting two birds with a stone. 143 144 145 146 147 ........................................................................... 148 149 150 151 4/3/97 rank and shape of a flag input. By default, the module builder was using the test 152 153 if (type != TYPE_INT || category != CATEGORY_REAL || rank != 1 || shape != 2) 154 DXSetError(ERROR_DATA_INVALID, "input \"preserve_volume\""); 155 156 the right test to use is 157 158 if (type != TYPE_INT || category != CATEGORY_REAL || !((rank == 0) || ((rank == 1)&&(shape == 1)))) 159 160 161 162 ........................................................................... 163 164 165 166 4/3/97 DXGetComponentValue(field, "positional_error") 167 168 the "positional_error" component was ignored by the simplification because 169 of a faulty test for the "dep" attribute 170 if was using 171 172 if (strcmp("positions", DXGetString((String)attr))) 173 174 instead of 175 176 if (!strcmp("positions", DXGetString((String)attr))) 177 178 which is the "c" translation of : are "positions" and DXGetString((String)attr) the same 179 180 the (!) will presumably get me many more times in the future. This is the advantage of 181 documenting the bugs. Hopefully I'll remember this one. 182 183 184 185 ........................................................................... 186 187 188 4/3/97 bug in the routine _dxfFloatDataBoundingBoxDiagonal. I noticed it when running 189 an example that fed several isosurfaces to the module. The bounding box value was inconsistent. 190 191 the code was as follows: 192 193 for(j=0;j<data_dim;j++) 194 { 195 bbox[j] = bbox[j+data_dim] = data[j]; 196 } 197 198 199 cur_data = data + data_dim; 200 201 for (i=0;i<n_data;i++,cur_data+=data_dim) 202 203 for(j=0;j<data_dim;j++) 204 { 205 if (cur_data[j] < bbox[j]) bbox[j] = cur_data[j]; 206 207 else if (cur_data[j] > bbox[j+data_dim]) bbox[j+data_dim] = cur_data[j]; 208 } 209 210 The problem was that the first set of data values was used to initialize the bounding box 211 and *then* after incrementing cur_data, I was examining n_data values. 212 Hence I was peeking outside of the data array and picking any ridiculous value there or 213 risking a core dump 214 215 one solution is to do the loop as 216 217 for (i=1;i<n_data;i++,cur_data+=data_dim) 218 219 220 ........................................................................... 221 222 223 4/3/97 GET_POSITIONAL_ERROR the existing positional error was not used. I noticed it because 224 the module was apparently simplifying too much on a second simplification. 225 226 what happened is that the portion of the code 227 228 starting with 229 230 array = (Array)DXGetComponentValue(field, "positional_error"); 231 232 was just before a goto label GET_OTHER_INPUTS, and the code jumped to that goto in case 233 no data "dep" positions was available, which concerns most of the cases. 234 235 so I changed the name of the goto label to GET_POSITIONAL_ERROR and placed it before the 236 relevant code. 237 238 ........................................................................... 239 240 241 4/3/97 this was an important omission rather than a bug. 242 243 when running the debugged module on non-manifold input, I realize that if 244 some faces in the input were degenerate, (if the three vertices of the face 245 weren't all different), the conversion to manifold would not eliminate such faces. 246 247 Curiously, this elimination of degenerate faces was done in the first version 248 of the conversion to manifold, but when re-programming it I decided that it 249 was not necessary. I think that it *is* necessary, because it perturbs 250 the joining of the corners. I am supposed to join corner c1 and c2 and c3 and c4 251 and all the corners must be different. If a face is degenerate, then for 252 instance c1 and c3 are the same and I don't know what will happen. 253 254 This is a nasty omission, because a lot of code was depending on the number of 255 faces not being changed in the conversion to manifold. Now there is a "face_lut" 256 and a lot of code needed to be updated. This took approximately 3 hours. 257 258 However, after incorporating this capability, problems still subsist, 259 leading us to the next bug: 260 261 262 ........................................................................... 263 264 265 4/3/97 - 4/4/97 dealing with erroneous input 266 267 in order to convert non manifold input surfaces to manifold surfaces, 268 it is necessary to identify the manifold edges and to join triangles 269 that share a manifold edge. If the edge is not manifold, no triangles can 270 be joined there. 271 272 ........................................................................... 273 274 4/4/97 with new method for fixing non-manifold input, the edge count 275 was not correct. I was counting too many edges and hence there were some 276 additional spurious edges. 277 278 In fact, in the routine _dxfLabelEdges() it is absolutely *necessary* 279 to compute the number of edges again, as some triangles might be 280 implicitely joined and there may be fewer edges than we previously thought. 281 282 This bug was very difficult to debug because the -g version was running fine 283 in the debugger, but I was convinced that it was a problem of memory being overwritten. 284 In practice, no memory was overwritten but there were fewer valid elements in the 285 edges array than the counter predicted leading in any catastrophic behavior when 286 such edges were processed. 287 288 ........................................................................... 289 290 4/4/97 SIGH! as if it wasn't enough for today: I fixed the counting problem for 291 the number of edges nE, and thus passed the pointer &nE to a bunch of routines 292 293 I forgot that _dxfInitHashTable(e_table, _dxfFlog2(nE) + 1,... 294 requires an estimate of the number of edges ! nE. 295 It was receiving a NE of 0!!!, resulting in all the hash entries being added to the 296 same list. 297 298 I noticed it because the program was soooo slooooow! 299 300 ........................................................................... 301 302 4/6/97 definitions of uint, u_int, ushort, u_short ... 303 304 the problem is that if a compiler knows about them on a particular 305 platform, another compiler on a different platform might not know! 306 307 below is a series of directives ifdef and include that seem to 308 work 309 310 ........................................................................... 311 312 4/7/97 DXArrayConvert again 313 314 I was getting an error ERROR: SimplifySurface: Invalid data: data ranks dont match 315 that after one hour of debugging, I attributed to DXArrayConvert again 316 so now I call instead of DXArrayConvert(array, TYPE_FLOAT, CATEGORY_REAL, 0) 317 318 data_array = DXArrayConvert(array, TYPE_FLOAT, CATEGORY_REAL, rank); 319 320 ........................................................................... 321 322 4/7/97 323 324 line 6226 of the _simplifysurface.c file 325 326 the line was 327 328 s_comp[i] = s_comp[star1[i1]]; 329 330 instead of: 331 332 s_comp[i] = t_comp[star1[i1]]; 333 334 resulting in reading the memory in some strange place 335 it took me 3 hours to hunt this bug which was causing a segmentation 336 fault on a dec machine (but not on other machines apparently) 337 the debugger dbx was very useful for that, as I could notice that 338 the array s_comp was not filled properly. It was surprisingly hard 339 to finally find this typo, though, who easily survived compiling. 340 341 ........................................................................... 342 343 4/7/97 344 345 *new_edges = (EdgeS *) DXReAllocate(*new_edges, nE * sizeof (EdgeS)); 346 line 760 of _simplifysurface.c 347 348 the problem is that I was reallocating *after* filling up the h-table 349 of edges. The edges were copied to another place and the edge entries 350 in the h-table were pointing to nothing=or to a completely different 351 piece of memory. 352 353 This bug only appeared in the purified version of the dxexec because 354 the reallocate was a realloc with a smaller size, which results in 355 leaving the memory where it is with most implementations, except 356 the implementation of Paula for the memorystubs.o which mallocs, 357 copies and frees (places the information in a completely different memory 358 location. 359 360 The most annoying part of this bug is that I am sure to have encountered 361 it before. I should have diagnosed it much sooner. 362 363 364 ........................................................................... 365 366 4/8/97 367 368 SQRT: DOMAIN ERROR observed on hp and on sun. What I did was to 369 change most of the sqrt calls, except those embedded into macros, 370 and make sure that do not call sqrt on 0 or negative input. 371 372 ........................................................................... 373 374 4/9/97 375 376 The "positional_error" component was probably resampled like 377 the other components. Then there should have been two components 378 called "positional_error" dep on "positions": the newly created component 379 + the resampled component. Apparently there was only one?? 380 381 What I did was to add a test to make sure that "positional_error" 382 would not be resampled. 383 384 ........................................................................... 385 386 4/10/97 387 388 function 389 Nodetype *_dxfnewnode() 390 391 I wasnt checking that what was returned was different from NULL 392 before doing 393 394 Nodetype *retval = _dxfgetnode(table); 395 396 retval->code = code; 397 retval->rec = rec; 398 399 now it reads 400 if (retval) 401 { 402 retval->code = code; 403 retval->rec = rec; 404 } 405 return(retval); 406 407 I hope that it will get rid of the "Memory Corrupt" signals 408 409 ........................................................................... 410 411 4/11/97 improvement rather than bug fix 412 413 I added a "data" switch in order to turn on or off the use of 414 vertex related data to decide whether to collapse edges 415 ........................................................................... 416 417 4/12/97-4/13/97 improvement rather than bug fix 418 419 I added a "stats" switch in order to provide information on 420 how many vertices and triangles the 421 module starts with and how many vertices and triangles the module ends up 422 with 423 424 425 ........................................................................... 426 427 428 4/15/97 Lloydt reported that after SimplifySurface, various attributes 429 of dx fields were lost 430 431 --> call DXCopyAttributes() after the components are created 432 433 when printing the objects resulting from SimplifySurface, I noticed 434 that sometimes the Bounding Box and Neighbors were created and sometimes 435 they werent. 436 437 --> call DXEndField() after the field is created (after the worker 438 routine has completed 439 440 441 442 ........................................................................... 443 444 445 4/30/97 Lloydt noted that SimplifySurface passes "quads" through 446 without processing them 447 448 --> implemented a _dxfQuads2Triangles internal conversion routine 449 from quads to triangles. The same connection dependent 450 data mapping can be used as before using the "face_lut" array. 451 452 453 454 ........................................................................... 455 456 4/30/97 treat the following classes in the doLeaf() handler 457 458 case CLASS_PRIVATE: 459 case CLASS_LIGHT: 460 case CLASS_CAMERA: 461 462 (before, I would have set an error "ERROR_BAD_CLASS", "encountered in object traversal"); 463 464 ........................................................................... 465 466 4/30/97 series are handled better 467 468 --> set is_series = (DXGetGroupClass((Group)in[0]) == CLASS_SERIES); 469 470 471 --> use 472 473 new_in[0] = DXGetSeriesMember((Series)in[0], i, &position); 474 475 DXSetSeriesMember((Series)new_group, i, (double) position, new_out[0]); 476 477 ........................................................................... 478 479 note: most of the 4/29/96 changes were done during the vacation period 480 4/15--4/28. They were checked on 4/30/96 481 482 ........................................................................... 483 484 5/1/97 485 486 I discovered that the "memory corrupt" signals are set by the DXFree routine 487 when it frees a portion of memory that was already freed. 488 489 in the routine _dxfBuildOrientedManifold I had the following piece of code: 490 491 if (*vlut) 492 {DXFree((Pointer)(*vlut));} 493 494 *vlut = new_vlut; 495 496 be careful, here, in case there is an error, new_vlut will be freed 497 498 499 I had to add the following after the error: flag 500 501 if (new_vlut) 502 { 503 DXFree((Pointer)new_vlut); 504 new_vlut was computed and vlut freed and replaced with new_vlut, 505 so we must set *vlut = NULL 506 *vlut = NULL;} 507 508 509 also I redid the routines _dxfCreateSimplificationDataStructure 510 _dxfFreeSimplificationDataStructure 511 _dxfPartialFreeSimplificationDataStructure to make sure that pointers were set 512 to NULL after being freed. 513 514 _dxfPartialFreeSimplificationDataStructure was created to make space and to avoid running 515 out of memory in _dxfBuildSimplifiedSurface 516 517 also. In 518 _dxfBuildSimplifiedSurface I made sure pointers were set to NULL after being 519 freed, to avoid freeing them twice in the error: part of the code 520 521 the building of new surface vertices and new surface triangles was separated 522 for clarity and for freeing un-used memory before allocating new memory 523 524 ........................................................................... 525 526 5/7/97 bug of Floating Point exception noticed on DEC 527 528 in the function _dxfSimplifiedBoundaryVertexLocation() 529 530 around lines 3700-3800 of _simplifysurface.c 531 532 A special case is encountered if c == 0, 533 then we also have x = 0 and the formula is simply ye=a 534 this case is handled with the general case without particular 535 attention, except that we can't divide by c and 536 also we should not execute the portion of code relative to 537 t == 0 or t == 1 538 539 ........................................................................... 540 */ 541 542 543 #include <stdio.h> 544 #include <math.h> 545 546 #if defined(HAVE_VALUES_H) 547 #include <values.h> 548 #endif 549 550 #if defined(HAVE_SYS_BSD_TYPES_H) 551 #include <sys/bsd_types.h> 552 #endif 553 554 #if defined(HAVE_SYS_TYPES_H) 555 #include <sys/types.h> 556 #endif 557 558 /* some vector operations: */ 559 560 #undef SQUARED_DISTANCE3 561 #define SQUARED_DISTANCE3(u,v) (((u)[0]-(v)[0])*((u)[0]-(v)[0]) + ((u)[1]-(v)[1])*((u)[1]-(v)[1]) +\ 562 ((u)[2]-(v)[2])*((u)[2]-(v)[2])) 563 564 565 /* At some point, I might wast to make sure that sqrt is called on non-zero input. 566 I will have to replace the macro with a function call */ 567 568 #undef DIST3 569 #define DIST3(u,v) ( sqrt ( SQUARED_DISTANCE3((u),(v)))) 570 571 #undef SCALPROD3 572 #define SCALPROD3(u,v) ( (u)[0]*(v)[0] + (u)[1]*(v)[1] + (u)[2]*(v)[2] ) 573 574 #undef NORM3 575 #define NORM3(u) ( sqrt ( SCALPROD3((u),(u)) ) ) 576 577 #undef MakeVec 578 #define MakeVec(u,v,w) {(w)[0] = (v)[0] - (u)[0]; (w)[1] = (v)[1] - (u)[1]; (w)[2] = (v)[2] - (u)[2];} 579 580 #undef AddVec 581 #define AddVec(u,v) {/* u = u + v */ (u)[0] += (v)[0]; (u)[1] += (v)[1]; (u)[2] += (v)[2];} 582 583 #undef VecCopy 584 #define VecCopy(u,v) {(v)[0] = (u)[0]; (v)[1] = (u)[1]; (v)[2] = (u)[2];} 585 586 #undef VecProd3 587 #define VecProd3(u,v,w) { (w)[0] = (u)[1]*(v)[2]-(u)[2]*(v)[1]; (w)[1] = (u)[2]*(v)[0]-(u)[0]*(v)[2]; \ 588 (w)[2] = (u)[0]*(v)[1]-(u)[1]*(v)[0];} 589 590 #undef EdgeBarycentricCoordinates 591 #define EdgeBarycentricCoordinates(u,v,a,b) { a = SCALPROD3(u,v); b = SCALPROD3(u,u); } 592 593 /* prototype of the driver routine for simplification */ 594 595 int _dxfSimplifySurfaceDriver(int nV, float *v, int nT, int *t, 596 int dim_data, 597 float *vertex_data, 598 float tolerance, 599 float data_tolerance, 600 /* to be able to resample data defined on a non-manifold 601 surface, we need look-up tables between vertices before 602 and after manifold conversion and between triangles before 603 and after manifold conversion */ 604 int *old_nV_after_manifold_conversion, 605 int *old_nT_after_manifold_conversion, 606 int *new_nV, /* new number of vertices after simplification*/ 607 float **new_v, /* new vertex array (positions) after 608 simplification */ 609 float **new_vertex_data, /* new data values after simplification */ 610 int *new_nT, /* new number of triangles after simplification*/ 611 int **new_t, /* new triangle array (connections) after 612 simplification */ 613 int **vertex_parents, /* relationship between simplified vertices 614 and manifold converted vertices */ 615 int **face_parents, /* relationship between simplified triangles 616 and manifold converted triangles */ 617 int **vertex_lut, /* for conversion from non-manifold */ 618 int **face_lut, /* to manifold surfaces */ 619 620 float *old_positional_error, 621 622 /* in case we resimplify a surface that was already simplified, 623 we use the old positional error as a starting point. The old positional 624 error may also have been provided by uncertainty measurements */ 625 626 float **new_positional_error, 627 float **face_normals, /* normals of the simplified faces */ 628 float **old_face_areas, /* face areas after conversion to manifold */ 629 float **new_face_areas, /* face areas after simplification */ 630 int preserve_volume, /* flag = 1 is the volume should be preserved */ 631 int simplify_boundary, /* flag = 1 if the boundary should be 632 simplified */ 633 /* flag = 1 if the boundary length should be preserved, in case 634 the boundary is simplified (simplify_boundary = 1) */ 635 int preserve_boundary_length); 636 637 638 639 float _dxfFloatDataBoundingBoxDiagonal(float *data, int n_data, int data_dim); 640 /* compute the bounding box of any multidimensional array of data */ 641 642 643 int _dxfVertexNormalsfromTriangleNormals(Array nor, int nV, int nT, int *triangles, 644 float *face_normals, float *face_areas); 645 /* convert the triangle normals to vertex normals */ 646 647 /* resample components dep positions and dep connections after simplification 648 "resample" means: 1) create a new component with the same name for the simplified surface 649 2) compute new value at each new vertex or at each new triangle 650 by averaging the values at the corresponding old vertices and the corresponding 651 old triangles: 652 653 the correspondences between old elements and new elements are obtained by 654 querying the following arrays: *vertex_parents, *face_parents, *vertex_lut, 655 *face_lut 656 657 the array "vertex_parents" gives for each old vertex after conversion to manifold 658 the new vertex to which it maps 659 the array "vertex_lut" gives for each vertex after conversion the corresponding 660 vertex before conversion. 661 662 the array "face_parents" gives for each old triangle after conversion to manifold 663 the new triangle to which it maps 664 the array "face_lut" gives for each triangle after conversion the corresponding 665 triangle before conversion. 666 667 When averaging vertex attributes, each vertex has weight 1 668 When averaging triangle attributes, each triangle has a weight equal to its area 669 670 the following components are not resampled 671 POSITIONS, CONNECTIONS, INVALID_POSITIONS INVALID_CONNECTIONS 672 NORMALS, NEIGHBORS, POSITIONAL_ERROR 673 and DATA if it is used by SimplifySurface. 674 675 */ 676 677 int _dxfResampleComponentValuesAfterSimplification(Field original_surface, 678 Field *simplified_surface, 679 int old_nV_after_conversion, 680 int old_nT_after_conversion, 681 int new_nV, 682 int new_nT, 683 int *vertex_parents, int *face_parents, 684 int *vertex_lut, int *face_lut, float *face_areas); 685 686 /* routine specialized for resampling position dependent attributes */ 687 688 Array _dxfResampleComponentDepPosAfterSimplification(Array array, 689 int old_nV_after_conversion, 690 int new_nV, int *vertex_parents, 691 int *vertex_lut); 692 693 /* routine specialized for resampling connection dependent attributes */ 694 695 Array _dxfResampleComponentDepConAfterSimplification(Array array, 696 int old_nT_after_conversion, 697 int new_nT, 698 int *face_parents, int *face_lut, 699 float *face_areas); 700 701 702 /*+--------------------------------------------------------------------------+ 703 | | 704 | _dxfComputeVertexValence | 705 | | 706 +--------------------------------------------------------------------------+*/ 707 708 /* using the set of triangles, compute the valence of each vertex */ 709 710 #define _dxfComputeVertexValence(nT, t, valence) \ 711 { \ 712 register int i, i1; \ 713 register u_char k; \ 714 \ 715 /* for each triangle, increment the valence of its three vertices */ \ 716 \ 717 for (i=i1=0; i<(nT); i++, i1+=3) for (k=0; k<3; k++) (valence)[(t)[i1+k]]++; \ 718 } 719 720 721 /* 722 following the suggestion of Greg Abram 723 724 I am using macros to implement the resampling functions, 725 so that they can be used on all the different scalar data types 726 without having to copy the code every time. 727 728 A few other macros are defined later on to reduce the number of 729 function calls on some computation-intensive routines 730 */ 731 732 733 734 /*+--------------------------------------------------------------------------+ 735 | | 736 | _dxfRESAMPLE_DEP_POS_AFTER_SIMPLIFICATION | 737 | | 738 +--------------------------------------------------------------------------+*/ 739 740 #define _dxfRESAMPLE_DEP_POS_AFTER_SIMPLIFICATION(type) \ 741 { \ 742 register int j; \ 743 \ 744 type \ 745 *old_data = (type *) DXGetArrayData(array), \ 746 *old_data_current, \ 747 *new_data = (type *) new_array_data; \ 748 \ 749 float \ 750 *vertex_data_averages_current; \ 751 \ 752 int \ 753 total_new_num_elements = new_nV * tensor_size; \ 754 \ 755 /*for each data element (including each tensor element) \ 756 we compute the average using the vertex lut and vertex parent \ 757 \ 758 we need a vertex_weights array allocated to zero beforehand \ 759 \ 760 we need an array of floats to store temporarily the averages \ 761 at each vertex \ 762 */ \ 763 \ 764 \ 765 for(i=0; i<old_nV_after_conversion; i++) /* this is the old number of vertices \ 766 after the conversion to manifolds */ \ 767 { \ 768 int \ 769 \ 770 old_vertex = vertex_lut[i], \ 771 new_vertex = vertex_parents[i]; \ 772 \ 773 /* we add one to the vertex weights of the parents */ \ 774 \ 775 vertex_weights[new_vertex] += 1; \ 776 \ 777 /* we add the tensor_size values to the average */ \ 778 \ 779 /* get the right offset in both old data and new data arrays */ \ 780 \ 781 old_data_current = old_data + old_vertex * tensor_size; \ 782 \ 783 vertex_data_averages_current = vertex_data_averages + new_vertex * tensor_size; \ 784 \ 785 for (j=0;j<tensor_size;j++) \ 786 \ 787 vertex_data_averages_current[j] += (float) old_data_current[j]; \ 788 \ 789 } \ 790 \ 791 /* then loop on all the averages, and divide them by the vertex_weights */ \ 792 \ 793 vertex_data_averages_current = vertex_data_averages; \ 794 \ 795 for(i=0; i< new_nV; i++,vertex_data_averages_current+=tensor_size) \ 796 \ 797 for (j=0;j<tensor_size;j++) \ 798 \ 799 vertex_data_averages_current[j] /= vertex_weights[i]; \ 800 \ 801 /* finally convert the vertex_data_averages to the requested scalar type */ \ 802 \ 803 \ 804 for(i=0; i< total_new_num_elements; i++) \ 805 \ 806 new_data[i] = (type) vertex_data_averages[i]; \ 807 } 808 809 810 /*+--------------------------------------------------------------------------+ 811 | | 812 | _dxfRESAMPLE_DEP_CON_AFTER_SIMPLIFICATION | 813 | | 814 +--------------------------------------------------------------------------+*/ 815 816 #define _dxfRESAMPLE_DEP_CON_AFTER_SIMPLIFICATION(type) \ 817 { \ 818 register int j; \ 819 \ 820 type \ 821 *old_data = (type *) DXGetArrayData(array), \ 822 *old_data_current = old_data, \ 823 *new_data = (type *) new_array_data; \ 824 \ 825 float \ 826 *face_data_averages_current; \ 827 \ 828 int \ 829 total_new_num_elements = new_nT * tensor_size; \ 830 \ 831 for(i=0; i< old_nT_after_conversion; i++, old_data_current += tensor_size) \ 832 /* this is the old number of faces after conversion and before simplification */ \ 833 { \ 834 int \ 835 old_face = face_lut[i], \ 836 new_face = face_parents[i]; \ 837 \ 838 /* we add the area of the face to the weight of the parent triangle */ \ 839 \ 840 face_data_weights[new_face] += old_face_areas[old_face]; \ 841 \ 842 /* we next update the averages for all associated tensor_size values */ \ 843 \ 844 /* get the right offset in both old data and new data arrays */ \ 845 \ 846 old_data_current = old_data + old_face * tensor_size; \ 847 \ 848 face_data_averages_current = face_data_averages + new_face * tensor_size; \ 849 \ 850 for (j=0;j<tensor_size;j++) \ 851 \ 852 face_data_averages_current[j] += old_face_areas[old_face] * (float) \ 853 old_data_current[j]; \ 854 \ 855 } \ 856 \ 857 /* then loop on all the averages, and divide them by the face_weights */ \ 858 \ 859 face_data_averages_current = face_data_averages; \ 860 \ 861 for(i=0; i< new_nT; i++,face_data_averages_current+=tensor_size) \ 862 \ 863 for (j=0;j<tensor_size;j++) \ 864 \ 865 face_data_averages_current[j] /= face_data_weights[i]; \ 866 \ 867 /* finally convert the vertex_data_averages to the requested scalar type */ \ 868 \ 869 \ 870 for(i=0; i< total_new_num_elements; i++) \ 871 \ 872 new_data[i] = (type) face_data_averages[i]; \ 873 } 874 875 876 /* copy an array from one type to another type */ 877 878 /*+--------------------------------------------------------------------------+ 879 | | 880 | CopyType2Type | 881 | | 882 +--------------------------------------------------------------------------+*/ 883 884 #define CopyType2Type(type1, type2, a, b, n) \ 885 { \ 886 register int i; \ 887 \ 888 type1 *a1 = (type1 *) (a); \ 889 \ 890 type2 *b1 = (type2 *) (b); \ 891 \ 892 for (i=0; i<(n) ; i++) b1[i] = (type2) a1[i];} 893 894 895 /* convert type "quad" connections to type "triangles" */ 896 897 int _dxfQuads2Triangles(int *nT, int **connections, int **face_lut); 898 899 900 /* data structures for the surface edges and associated routines */ 901 902 903 typedef struct _EdgeS_ { 904 905 short label; /* used for a flag or for 906 classifying edges among boundary, regular and singular edges */ 907 int v[2]; /* end points */ 908 int t[2]; /* incident triangles */ 909 } EdgeS; 910 911 /* 912 surface edges have an orientation: 913 v[0] >= v[1] and t[0] is before t[1] in clockwise order 914 when we rotate around v0. 915 916 /|\v1 917 / | \ 918 / | \ 919 / t0|t1 \ 920 \ | / 921 \ | / 922 \ | / 923 \|/v0 v0 >= v1 924 */ 925 926 /* routines for indexing edges in a hash table */ 927 928 int _dxfEdgeSKey(EdgeS *); 929 int _dxfEdgeSFilter(EdgeS *, EdgeS *); 930 931 932 #define SurfEdgeLabel(edg) (edg)->label 933 #define SurfEdgeGetVertices(edg) (edg)->v 934 #define SurfEdgeGetTriangles(edg) (edg)->t 935 936 937 /* my hash tables: 938 939 they do not have the restriction that no more than 940 16 elements should have the same code, that is documented in the DX hash tables: 941 942 if someone gives me a surface with 1,000,000 triangles, it should 943 roughly have 1,500,000 edges and I don't know of a simple coding 944 function that would never give 16 times the same code to edges (it 945 may be possible to design such a function but I don't know how) 946 947 */ 948 949 /*...........................................................................*/ 950 951 /* Following are includes needed for hash-tables using linear probing. 952 In case of a code conflict, data items are chained */ 953 954 /*...........................................................................*/ 955 956 typedef char *RECTYPE ; 957 typedef char *KEYTYPE ; 958 959 typedef struct _nodetype_{ 960 int code; 961 RECTYPE rec; 962 struct _nodetype_ *next; 963 } Nodetype; 964 965 typedef struct _Block_{ 966 RECTYPE rec; 967 struct _Block_ *next; 968 } Block; 969 970 typedef struct _htable_{ 971 Nodetype **bucket; 972 Block *block; 973 Nodetype *freenode; 974 975 int nodes_per_block ; 976 int size, mask; 977 int search_code; 978 Nodetype *search_node; 979 KEYTYPE search_key; 980 981 int (*fcode)(KEYTYPE); 982 int (*filter)(RECTYPE, RECTYPE); 983 KEYTYPE (*compute_key)(RECTYPE); 984 985 } Htable; 986 987 /*...........................................................................*/ 988 989 /* prototypes for H-table procedures 990 */ 991 992 /*...........................................................................*/ 993 994 /* add a record to the hash-table */ 995 int _dxfAddRecord(Htable *, RECTYPE, KEYTYPE); 996 KEYTYPE _dxfIdentifyRecord2Key(RECTYPE rec); 997 998 999 /* initialize the hash-table */ 1000 1001 int _dxfInitHashTable(Htable *, int, int (*fcode)(KEYTYPE), 1002 int (*filter)(RECTYPE,RECTYPE), 1003 KEYTYPE (*compute_key)(RECTYPE)); 1004 1005 /* determine whether a particular record is in the hash-table and if yes, find that record 1006 meaning, return a pointer to that record*/ 1007 1008 char *_dxfFindRecord(Htable *, KEYTYPE); 1009 void _dxfHashTableReset(Htable *); 1010 1011 int _dxfNewBlock(Htable *); 1012 1013 /* once _dxfFindRecord() has been called, retrieve the next record with the same code and 1014 passing through the same filter() from the hash table */ 1015 1016 char *_dxfNextRecord(Htable *); 1017 Nodetype *_dxfgetnode(Htable *); 1018 Nodetype *_dxfnewnode(Htable *,RECTYPE, int); 1019 1020 int _dxfFlog2(int n); /* log base two using integer arithmetic */ 1021 1022 /* reset the hash table and free the buckets 1023 after that, either (1) inithashtable must be called again to refill the hash table with new records or 1024 completely different data 1025 1026 or (2) free(table) to finish it off 1027 unless it was defined as Htable table[1] in which case nothing more has to be done */ 1028 1029 #define _dxfFreeHashTable(table) {_dxfHashTableReset(table);if ((table)->bucket) \ 1030 {DXFree((Pointer)(table)->bucket);(table)->bucket = NULL;}} 1031 1032 1033 /* find a particular surface edge in a Hash table indexing edges */ 1034 1035 EdgeS *_dxfFindEdge(Htable *e, int vA, int vB); 1036 1037 1038 1039 /* convert a surface mesh to an oriented manifold surface mesh 1040 manifold means that each vertex has a single sheet of triangles touching it 1041 and each edge has at most two incident triangles */ 1042 1043 int _dxfToOrientedManifold(int nV, float *v, int nT, int *t, 1044 int *nVm, float **vm, int *nTm, int **tm, int **vlut, int **flut, 1045 int *n_edges, Htable *e_table, EdgeS **e); 1046 1047 /* this is the first operation of conversion to a manifold: 1048 triangles that are degenerate are eliminated 1049 a triangle is degenerate if either (1) one or more vertices are less then 0 or more than nV 1050 or (2) two or three vertices in the triangle are the same 1051 */ 1052 1053 int _dxfEliminateDegenerateTriangles(int nV, int nT, int *t, int *new_nT, int **new_t, int **flut); 1054 1055 /* second operation of conversion to manifold: 1056 the vertices that are not referenced after eliminating the degenerate triangles are 1057 eliminated from the list of vertices */ 1058 1059 int _dxfEliminateStandaloneVertices(int nV, float *v, int *nVm, float **vm, int nT, int *t, 1060 int **vlut); 1061 1062 /* then build a list of edges that are manifold edges, and join the corners of the triangles 1063 that share such an edge */ 1064 1065 int _dxfJoinTriangleCorners(int nT, int *t, Htable *e_table, int *n_edges, EdgeS **edges, 1066 int *fathers); 1067 1068 int _dxfBuildOrientedManifold(int *nV, float **v, int nT, int *t, int **vlut, 1069 int *nE, EdgeS **edges, Htable *e_table, int *fathers); 1070 1071 int _dxfIndexEdges(int *nE, EdgeS **edges, Htable *e_table, int nT, int *t); 1072 1073 void _dxfInitFather(int n, int *father); 1074 1075 1076 /* perform path conpression on the representatives of a forest data structure 1077 the macro FATHER is more efficient than the _dxfFather(i,f) procedure 1078 it calls the procedure only if necessary */ 1079 1080 #define FATHER(i, f) (((i) == (f)[i]) ? (i) : \ 1081 ((f)[i] == (f)[(f)[i]]) ? (f)[i] : _dxfFather(i,f)) 1082 1083 int _dxfFather(int i, int *father); 1084 1085 int _dxfJoin(int i, int j, int *father); 1086 1087 1088 1089 /*...........................................................................*/ 1090 /* 1091 default parameters for simplification 1092 */ 1093 /*...........................................................................*/ 1094 1095 1096 #define SIMPLIFY_MAX_ANGLE_NOR 1.5 /* maximum angle between normals of corresponding triangles 1097 before and after an edge-collapse */ 1098 #define SIMPLIFY_COMPACTNESS_RATIO .2 /* maximum ratio between the minimum compactness 1099 in the vertex star after simplification and the 1100 minimum compactness in the edge star before simplification */ 1101 #define SIMPLIFY_VALENCE_MAX 30 /* maximum valence at a vertex after simplification */ 1102 #define SIMPLIFY_VALENCE_MAX1 31 1103 #define SIMPLIFY_VALENCE_MAX4 34 1104 1105 #define FOUR_SQRT_3 6.928203230275508 1106 #define EPS_VERTEX_OUTSIDE_TRIANGLE -1e5 /* threshold on the barycentric coordinate such that if 1107 a barycentric coordinate is less than that the corresponding 1108 point is decided to be on the other side of the edge */ 1109 1110 /* typedef for some simple geometric entities: */ 1111 1112 typedef float Vertex[3]; 1113 typedef float Plane[4]; 1114 typedef double VertexD[3]; 1115 typedef int Face[3]; 1116 1117 #define DIRECTION_COUNTER_CLOCKWISE 2 1118 #define DIRECTION_CLOCKWISE 1 1119 1120 /*...........................................................................*/ 1121 1122 /* 1123 data structures for binary heaps 1124 */ 1125 1126 /*...........................................................................*/ 1127 1128 1129 typedef struct _Heap_ { 1130 int n,last; 1131 int *perm; 1132 int *invp; 1133 float *key; 1134 1135 } Heap; 1136 1137 #define LeftSon(i) 2*(i)+1 1138 #define RightSon(i) 2*(i)+2 1139 1140 Heap *_dxfNewHeap(int n); 1141 unsigned long _dxfHeapSize(int n); 1142 void _dxfResetHeap(Heap *heap); 1143 int _dxfHeapAdd(float k, Heap *heap); 1144 int _dxfHeapDelMin(Heap *heap); 1145 int _dxfHeapDelete(int i, Heap *heap); 1146 1147 #define HeapLength(heap) (heap)->last 1148 #define HeapLengthMax(heap) (heap)->n 1149 #define HeapFull(heap) (HeapLength(heap) == HeapLengthMax(heap)) 1150 #define HeapValidIndex(heap, index) \ 1151 ((index) < HeapLengthMax(heap) && (index) >= 0) 1152 1153 /* you can't remove these elements from the heap but you can add them in: */ 1154 1155 #define HeapOutsideIndex(heap) HeapLengthMax(heap) + 1 1156 #define HeapInitialIndex(heap) HeapLengthMax(heap) 1157 1158 /* you can neither add nor remove these elements */ 1159 1160 #define HeapInvalidIndex(heap) -1 1161 1162 1163 /* internal structure used for simplification: */ 1164 1165 typedef struct _SimpData_ { 1166 1167 /* data marked "to be allocated" is managed by the 1168 1169 _dxfCreateSimplification.... 1170 1171 _dxfFreeSimplification.... 1172 1173 routines. it is used only internally to the _dxfSimplifySurfaceDriver routine and is not exported 1174 1175 */ 1176 1177 Heap 1178 *edge_heap; /* to be allocated */ 1179 1180 int 1181 nV, nT, nE, 1182 data_dim, 1183 *triangles, 1184 *edge2index, /* to be allocated */ 1185 *index2edge, /* to be allocated */ 1186 *vertex_father, 1187 *edge_father, /* to be allocated */ 1188 *triangle_father, 1189 *vertex_lut, 1190 preserve_volume, 1191 simplify_boundary, 1192 preserve_boundary_length, 1193 valence_max, 1194 1195 /* various counters: */ 1196 1197 num_edg_remaining_4_testing, 1198 num_edg_collapsed, 1199 num_edg_weights, 1200 edg_rejected_4_topology, 1201 edg_rejected_4_geometry, 1202 edg_rejected_4_tolerance, 1203 last_edge_added2heap, /* edge number for the last edge added to the heap */ 1204 heap_initial_index, 1205 heap_outside_index; 1206 1207 EdgeS 1208 *edge_array; 1209 1210 Htable 1211 *e_table; 1212 1213 Vertex 1214 *vert, 1215 *normal, /* to be allocated */ 1216 simplified_vertex; 1217 1218 u_char 1219 *boundary_vert; /* to be allocated */ 1220 1221 u_short 1222 *valence; /* to be allocated */ 1223 1224 float 1225 tolerance, 1226 data_tolerance, 1227 *err_volume, /* to be allocated */ 1228 *tol_volume, /* to be allocated */ 1229 *old_face_areas, 1230 *area, /* to be allocated */ 1231 *compactness, /* to be allocated */ 1232 min_scalprod_nor, 1233 compactness_ratio, 1234 *vx_data, /* to be allocated */ 1235 *vx_data_error, /* to be allocated */ 1236 *vx_data_potential_values, /* to be allocated */ 1237 /* potential data values at the simplified vertex */ 1238 *vx_old_data, 1239 *vx_new_data, /* necessary to compute the discrepancy in data values 1240 at the vertices (or corners) of the mapping between 1241 original and simplified surfaces */ 1242 1243 vx_data_potential_err; /* potential data error at the simplified vertex */ 1244 1245 int (*get_lowest_weight_edge)(struct _SimpData_ *); /* pointer to a function used 1246 for extracting the edge 1247 with lowest weight from 1248 the heap */ 1249 1250 int (*add_edge) (struct _SimpData_ *, int); 1251 int (*remove_edge) (struct _SimpData_ *, int, int); 1252 1253 /* other buffers used for storing information relative to vertex and edge stars */ 1254 1255 int 1256 *vertex_star_buffer; 1257 1258 Pointer 1259 edge_star_buffer_fl, 1260 edge_star_buffer_vx; 1261 1262 } SimpData; 1263 1264 1265 1266 int _dxfRemoveLowestWeightEdgeFromHeap(SimpData *simp_data); 1267 1268 int _dxfAddEdge2Heap(SimpData *simp_data, int the_edge); 1269 1270 int _dxfRemoveEdgeFromHeap(SimpData *simp_data, int index, int the_edge); 1271 1272 int _dxfBuildSimplifiedSurface(SimpData *simp_data, int *new_nV, float **new_v, 1273 int *new_nT, int **new_t, 1274 float **new_face_areas, float **face_normals, 1275 float **new_positional_error, float **new_vertex_data); 1276 1277 int _dxfSimplifyManifoldSurface(SimpData *simp_data); 1278 1279 int _dxfCreateSimplificationDataStructure(SimpData *simp_data, float *data, float *old_pos_err); 1280 1281 /* free all the simplification data structures: */ 1282 1283 int _dxfFreeSimplificationDataStructure(SimpData *simp_data, float *data, float *old_pos_err); 1284 1285 /* free only the simplification data structures that are not exported by the driver routine, 1286 to avoid running out of memory in _dxfBuildSimplifiedSurface() */ 1287 1288 int _dxfPartialFreeSimplificationDataStructure(SimpData *simp_data); 1289 1290 int _dxfTrianglesNormalAreaCompactness( int nT, Face *t, Vertex *v, Vertex *t_normal, 1291 float *t_area, float *t_compactness); 1292 1293 int _dxfTriangleNormalQR2D(VertexD *tri, double *n); 1294 1295 int _dxfTriangleNormalQR2(Vertex *tri, Vertex n); 1296 1297 int _dxfVectorProductQRD(VertexD x1, VertexD x2, double *n); 1298 1299 int _dxfFlagBoundaryVertices(int nE, EdgeS *edges, u_char *boundary_vert); 1300 1301 int _dxfMarkEdgesAdjacent2Boundary(SimpData *simp_data); 1302 1303 int _dxfBuildEdgeHeap(SimpData *simp_data); 1304 1305 void _dxfInitArray(char *array, char *data, int n, int size); 1306 1307 int _dxfCollapsibilityTestsBoundaryEdge2ndKind(SimpData *simp_data, 1308 int edge_num, int v0, int v1, int val0, int val1, 1309 int move_simplified_vertex); 1310 1311 int _dxfDirectedParentVertexStar(int parent_vertex, int first_triangle, 1312 int valence, int *star, SimpData *simp_data, int direction); 1313 1314 int _dxfRotateParentTriangle(int parent_vertex, int *vertex_fathers, 1315 int *tri_v, int *tri_v_parents, int *tri_v_rotated); 1316 1317 int _dxfRotateParentEdge(int parent_vertex, int *vertex_fathers, 1318 int *edge_v, int *edge_t, int *edge_v_parents, 1319 int *edge_v_rotated, int *edge_t_rotated); 1320 1321 /* macros used for finding the next edge when rotating around a vertex in 1322 on a simplified surface */ 1323 1324 1325 1326 #define _dxfNextParentEdge(oriented_tri, directionCW) \ 1327 _dxfFather( \ 1328 (int) (((EdgeS *)_dxfFindEdge(simp_data->e_table,(oriented_tri)[0],(oriented_tri)[directionCW]))\ 1329 -simp_data->edge_array), simp_data->edge_father) 1330 1331 1332 1333 #define _dxfNextParentTriangle(oriented_edge_t, directionCW) \ 1334 ((directionCW) == 1)? FATHER((oriented_edge_t)[1], simp_data->triangle_father):\ 1335 FATHER((oriented_edge_t)[0], simp_data->triangle_father) 1336 1337 int _dxfManifoldLinkTest(int *link0, int *link1, int val0, int val1); 1338 1339 int _dxfCmpIntSmallFirst(char *s1, char *s2); 1340 1341 float _dxfClosestPointOnEdge(Vertex X, Vertex XP, Vertex A, Vertex B, 1342 float *t);/* barycentric coordinates: t, 1-t */ 1343 1344 void _dxfMakeEdgeBarycenter(Vertex v0, Vertex v1, float alpha0, Vertex v); 1345 1346 float _dxfNormalize3Vector(Vertex v); 1347 1348 int _dxfSolveQuadraticEquation(double A, double B, double C, 1349 double *sol1, double *sol2, double eps); 1350 1351 int _dxfBoundaryCollapse2KndGeometricallyOk(SimpData *simp_data, int v0, int v1, 1352 int *star0, int *vstar0, int val0, Vertex *s_normal0, float *s_area0, 1353 int direction, float min_scalprod_nor, float compactness_ratio); 1354 1355 float _dxfFastCompactness3DTriangle2(Vertex *tri, float *triangle_area); 1356 1357 1358 1359 /*+--------------------------------------------------------------------------+ 1360 | | 1361 | _dxfTriangleBasisVectors | 1362 | | 1363 +--------------------------------------------------------------------------+*/ 1364 1365 1366 1367 #define _dxfTriangleBasisVectors(tri, x1, x2, origin) \ 1368 { \ 1369 /* compute two basis vectors of a triangle, such that the first basis vector \ 1370 will correspond to the shortest edge, \ 1371 return the number of the vertex that is chosen as the triangle origin */ \ 1372 \ 1373 \ 1374 int the_origin = 0; \ 1375 \ 1376 /* I- find the shortest edge of the triangle \ 1377 2 \ 1378 /\ \ 1379 edg1 / \ edg0 \ 1380 /____\ \ 1381 0 edg2 1 */ \ 1382 \ 1383 register u_char j; \ 1384 double min_len = SQUARED_DISTANCE3((tri)[0],(tri)[2]); \ 1385 u_char shortest_edge = 1; \ 1386 \ 1387 for (j=0;j<2;j++) \ 1388 { \ 1389 double len = SQUARED_DISTANCE3((tri)[j],(tri)[j+1]); \ 1390 \ 1391 if (len < min_len) \ 1392 { \ 1393 min_len = len; \ 1394 shortest_edge = (j+2)%3; \ 1395 } \ 1396 } \ 1397 \ 1398 the_origin = (shortest_edge + 1) %3; \ 1399 \ 1400 (origin) = the_origin; \ 1401 \ 1402 /* II- assign the smallest edge to one of the vectors x1 and x2 */ \ 1403 { \ 1404 u_char \ 1405 dest1 = (shortest_edge+2)%3, \ 1406 dest2 = shortest_edge; \ 1407 \ 1408 register u_char k; \ 1409 \ 1410 for(k=0;k<3;k++) \ 1411 { \ 1412 (x1)[k] = (tri)[dest1][k] - (tri)[the_origin][k]; \ 1413 (x2)[k] = (tri)[dest2][k] - (tri)[the_origin][k]; \ 1414 } \ 1415 } \ 1416 } 1417 1418 int _dxfErrorWithinToleranceVBoundary2ndKnd(SimpData *simp_data, int v0, int v1, 1419 int *vstar0, int *vstar1, int val0, int val1); 1420 1421 1422 #define ExtractTriangleFromVertexStar( vstar, i) a_triangle; { \ 1423 memcpy(a_triangle[1], simp_data->vert[(vstar)[(i)]], sizeof(Vertex)); \ 1424 memcpy(a_triangle[2], simp_data->vert[(vstar)[(i)+1]], sizeof(Vertex));} 1425 1426 /* macro for taking the minimum of the potential errors at the pivot vertex, 1427 also called the *simplified_vertex in my article */ 1428 1429 1430 1431 #define UPDATE_ERR_PIVOT \ 1432 { \ 1433 if (potential_err_pivot > largest_err_pivot) \ 1434 largest_err_pivot = potential_err_pivot; \ 1435 edge_collapsible = (largest_err_pivot <= tol_pivot);} 1436 1437 1438 1439 #define _dxfErrorTestDataOnSurfaceInit(data) \ 1440 { \ 1441 /* initialize the potential error value for the data at the potential simplified vertex location */ \ 1442 if ((data)->vx_data) {(data)->vx_data_potential_err = 0.0;}} 1443 1444 1445 1446 #define _dxfErrorTestDataOnSurfaceEnd(data,v) \ 1447 { \ 1448 \ 1449 /* replace the data error value at the simplified vertex with the potential data error value */ \ 1450 \ 1451 if ((data)->vx_data) \ 1452 {(data)->vx_data_error[(v)] = (data)->vx_data_potential_err; \ 1453 \ 1454 /* replace the data value with the potential data value */ \ 1455 memcpy((data)->vx_data + (v) * (data)->data_dim, (data)->vx_data_potential_values, \ 1456 (data)->data_dim * sizeof (float)); \ 1457 }} 1458 1459 1460 #define _dxfErrorTestDataOnSurfaceInterpolate2(data, v0, v1, t) \ 1461 { \ 1462 if ((data)->vx_data) \ 1463 { \ 1464 /* interpolate the data value at the potential simplified vertex */ \ 1465 \ 1466 int coordinate = 0; \ 1467 \ 1468 float \ 1469 *data0 = (data)->vx_data + (v0) * (data)->data_dim, \ 1470 *data1 = (data)->vx_data + (v1) * (data)->data_dim, \ 1471 t1 = 1. -(t); \ 1472 \ 1473 for (coordinate = 0; coordinate < (data)->data_dim; coordinate++) \ 1474 \ 1475 (data)->vx_data_potential_values[coordinate] = data0[coordinate] * (t) + data1[coordinate] * t1;\ 1476 } \ 1477 } 1478 1479 1480 #define _dxfErrorTestDataOnSurfaceInterpolate3(data, v0, v1, v2, t0, t1, t2) \ 1481 { \ 1482 if ((data)->vx_data) \ 1483 { \ 1484 /* interpolate the data value at the potential simplified vertex */ \ 1485 \ 1486 int coordinate = 0; \ 1487 \ 1488 float \ 1489 *data0 = (data)->vx_data + (v0) * (data)->data_dim, \ 1490 *data1 = (data)->vx_data + (v1) * (data)->data_dim, \ 1491 *data2 = (data)->vx_data + (v2) * (data)->data_dim; \ 1492 \ 1493 for (coordinate = 0; coordinate < (data)->data_dim; coordinate++) \ 1494 \ 1495 (data)->vx_data_potential_values[coordinate] = \ 1496 data0[coordinate] * (t0) + data1[coordinate] * (t1) + data2[coordinate] * (t2); \ 1497 } \ 1498 } 1499 1500 int _dxfErrorTestDataOnSurfaceUpdate(SimpData *simp_data, int new_v, int old_v, 1501 int new_v2, int new_v3, 1502 float alpha_sv, float alpha_n_v2, float alpha_n_v3, 1503 int old_v1, int old_v2, int old_v3, 1504 float alpha_o_v1, float alpha_o_v2, float alpha_o_v3); 1505 1506 float _dxfSegmentIntersection(Vertex A, Vertex B, Vertex C, Vertex D, 1507 float *lambda, float *mu, 1508 int changed_AorB, int changed_C, int changed_D); 1509 1510 int _dxfHouseholderPreMultiplication(float *A, int mA, float *v, int m, int n, float *w ); 1511 1512 int _dxfHouseholderPreMultiplicationDbl(double *A, int mA, double *v, int m, int n, double *w ); 1513 1514 int _dxfCollapseBoundaryEdge2ndKnd(SimpData *simp_data, int edg_num, int v0, int v1, int val0, 1515 int val1, int *star0, int *star1, int *vstar0, int *vstar1, 1516 int *estar0, int *estar1, Vertex *s_normal0, Vertex *s_normal1, 1517 float *s_area0, float *s_area1); 1518 1519 int _dxfRemoveEdgesFromHeap(int *estar, u_short val, int *edge2index, SimpData *simp_data); 1520 1521 int _dxfReinstateEdgesInHeap(int *estar, u_short val, SimpData *simp_data, 1522 int *num_edg_added); 1523 1524 int _dxfCollapsibilityTestsBoundaryEdge1stKind(SimpData *simp_data, 1525 int edge_num, int v0, int v1, int val0, int val1); 1526 1527 int _dxfManifoldLinkTest1stKnd(int v0, int val0, int *star1, int *link1, int val1, 1528 SimpData *simp_data); 1529 1530 int _dxfBoundaryCollapse1stKndGeometricallyOk(SimpData *simp_data, int v0, int *star0, int *vstar0, 1531 int val0, Vertex *s_normal0, float *s_area0, 1532 float min_scalprod_nor, float compactness_ratio); 1533 1534 int _dxfErrorWithinToleranceVBoundary1stKnd(SimpData *simp_data, int v0, int v1, int *vstar1, 1535 int val1); 1536 1537 int _dxfClosestPointOnTriangle(Vertex *tri, Vertex w, Vertex wp, Vertex bary, float *dist); 1538 1539 int _dxfTriangle3DBarycentricCoordinates2(Vertex *tri, Vertex w, Vertex wp, Vertex bary, 1540 float *residual); 1541 1542 void _dxfMakeBarycenter(int nV, Vertex *vv, Vertex bary, float *bary_coord); 1543 1544 int _dxfSolveSymmetric2x2Eqn(double a, double b, double d, double e, double f, 1545 double *x, double *y); 1546 1547 1548 int _dxfCollapseBoundaryEdge1stKnd(SimpData *simp_data, int edg_num, int v0, int v1, int val0, 1549 int val1, int *star1, int *vstar1, int *estar1, 1550 Vertex *s_normal1, float *s_area1); 1551 1552 int _dxfCollapseTopologicallyFeasible(SimpData *simp_data, int *vstar0, int *vstar1, 1553 u_short val0, u_short val1); 1554 1555 void _dxfBuildParentEdgeStars(EdgeS *edg, int v0, int v1, u_short val0, 1556 u_short val1, int *star0, int *star1, SimpData *s); 1557 1558 void _dxfParentVertexStar(int vf, int t0, u_short val, int *star, SimpData *simp_data); 1559 1560 void _dxfApplyTranslation(int nV, Vertex *vv, Vertex T); 1561 1562 void _dxfOppositeVector(Vertex u, Vertex v); 1563 1564 void _dxfCopyFaceNormalsAreasCompactness(Plane *plane, float *s_area, float *s_comp, 1565 int *star0, int *star1, int val0, int val013, 1566 Vertex *t_normal, float *t_area, float *t_comp); 1567 1568 float _dxfMinFloatArray(float *array, int n); 1569 1570 void _dxfSimplifiedVertexLocation(SimpData *simpdata, u_short val0, u_short val1, u_short star_val, 1571 Vertex *v, Face *f, Plane *p, float *area, float *w, 1572 u_char method); 1573 1574 double _dxfStarPlaneEquationsAndVolume(Vertex *star_vert, Face *star_face, Plane *star_plane, 1575 float *star_areas, int star_val); 1576 1577 void _dxfComposeVectors(Vertex u, Vertex v, float lambda, Vertex w); 1578 1579 int _dxfPositionVertexLSConstraints2(Plane *star_plane, u_short star_val, Vertex simplified_vertex); 1580 1581 void _dxfFPlaneForIdenticalVolume(u_short val0, u_short val1, Plane p, Vertex *star_vert, 1582 Face *star_face, float star_volume); 1583 1584 int _dxfCollapseGeometricallyFeasible(Vertex *svert, Face *sface, 1585 Plane *splane, float *old_comp, float new_comp_min, 1586 Vertex *snor0, Vertex *snor1, 1587 u_short val0, u_short sval, 1588 float min_scalprod, float compactness_ratio); 1589 1590 int _dxfNoFlipOverCheck(u_short first_face, u_short last_face, 1591 Vertex *nor, Plane *plane, float min_scalprod_nor); 1592 1593 int _dxfErrorWithinToleranceV(SimpData *simp_data, Vertex *star_vert, Face *star_face, 1594 Plane *star_plane, int *star0, int *star1, int val0, int val1); 1595 1596 void _dxfSimplifiedStarNormals(Face *sface, Vertex *svert, Vertex *snor, float *sarea, float *scomp, 1597 Vertex simpvert, u_short first_face, u_short last_face); 1598 1599 int _dxfCollapseEdge(SimpData *simp_data, int v0, int v1, int *star0, int *star1, 1600 u_short val0, u_short val1, Vertex *nor0, Vertex *nor1, float *area0, 1601 float *area1, float *comp0, float *comp1); 1602 1603 void _dxfReplaceFaceNormals(int *star, Vertex *new_nor, float *new_area, float *new_comp, 1604 Vertex *nor, float *area, float *comp, u_short val); 1605 1606 #endif /* _SIMPLESURF_H_ */ 1607 1608