/***********************************************************************/ /* Open Visualization Data Explorer */ /* (C) Copyright IBM Corp. 1989,1999 */ /* ALL RIGHTS RESERVED */ /* This code licensed under the */ /* "IBM PUBLIC LICENSE - Open Visualization Data Explorer" */ /***********************************************************************/ #include /* SimplifySurface is a module that approximates a triangulated surface and resamples data */ /*********************************************************************/ /* Author: Andre GUEZIEC, March/April 1997 */ /*********************************************************************/ #include "dx/dx.h" #include "simplesurf.h" /* this file contains the "guts" of the surface simplification procedure As much as possible I am using simple types For best memory management, I have converted my routines so that they use the data explorer DXAllocate and DXFree calls */ /*+--------------------------------------------------------------------------+ | | | _dxfSimplifySurfaceDriver | | | +--------------------------------------------------------------------------+*/ int _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, /* for fixing non-manifolds */ int **face_lut, /* same thing: for eliminating degenerate faces */ 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) /* all of the arrays that are ** must be allocated by the program if any of such arrays is not allocated properly, than the routine must return ERROR (zero). However, in case of an error, there is no need to DXFree() such arrays, as they are freed by the caller routine */ { int success = 0; int nE; EdgeS *e = NULL; Htable edge_table[1]; /* create a simplification data structure where all the relevant arrays are copied */ SimpData simp_data[1]; bzero(simp_data, sizeof (SimpData)); bzero(edge_table, sizeof (Htable)); /* the first step is to convert the surface to a manifold, in case Isosurface would have created non-manifold vertices (I know that happens sometimes) It is also possible that a surface provided directly by the user would not be a manifold */ if (_dxfToOrientedManifold(nV, v, nT, t, new_nV, new_v, new_nT, new_t, vertex_lut, face_lut, &nE, edge_table, &e)) { *old_nV_after_manifold_conversion = *new_nV; *old_nT_after_manifold_conversion = *new_nT; simp_data->nT = *new_nT; /* the number of vertices, edges and triangles */ simp_data->nE = nE; /* after conversion to manifold and before simplification */ simp_data->nV = *new_nV; simp_data->data_dim = (vertex_data) ? dim_data: 0; simp_data->preserve_volume = preserve_volume; simp_data->simplify_boundary = simplify_boundary; simp_data->preserve_boundary_length = preserve_boundary_length; simp_data->tolerance = tolerance; simp_data->data_tolerance = data_tolerance; *old_face_areas = (float *) DXAllocateZero( simp_data->nT * sizeof (float)); *vertex_parents = (int *) DXAllocateZero( simp_data->nV * sizeof (int)); *face_parents = (int *) DXAllocateZero( simp_data->nT * sizeof (int)); /* part of the simplification data structure that references existing arrays */ simp_data->edge_array = e; simp_data->e_table = edge_table; simp_data->vert = (Vertex *) *new_v; simp_data->triangles = *new_t; simp_data->vertex_father = *vertex_parents; simp_data->triangle_father = *face_parents; simp_data->old_face_areas = *old_face_areas; simp_data->vertex_lut = *vertex_lut; if (*old_face_areas && *vertex_parents && *face_parents && (_dxfCreateSimplificationDataStructure(simp_data, vertex_data, old_positional_error)) && (_dxfSimplifyManifoldSurface(simp_data))) { /* free only the simplification data structures that are not exported by the driver routine, to avoid running out of memory in _dxfBuildSimplifiedSurface(): */ _dxfPartialFreeSimplificationDataStructure(simp_data); /* build a simplified surface using the vertex_father and triangle_father arrays */ if (_dxfBuildSimplifiedSurface(simp_data, new_nV, new_v, new_nT, new_t, new_face_areas, face_normals, new_positional_error, new_vertex_data)) success = 1; _dxfFreeSimplificationDataStructure(simp_data, vertex_data, old_positional_error); } _dxfFreeHashTable(edge_table); if (e) DXFree((Pointer)e); } /* if the conversion to oriented manifold fails, vertex_lut, new_v, new_t ... etc. are freed in the caller routine */ return success; } /*+--------------------------------------------------------------------------+ | | | _dxfToOrientedManifold | | | +--------------------------------------------------------------------------+*/ int _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) { int success = 0; /* we first eliminate the degenerate triangles */ if (_dxfEliminateDegenerateTriangles(nV, nT, t, new_nT, new_t, flut)) { /* after that the number of triangles does not change any more */ nT = *new_nT; /* we then eliminate the standalone vertices */ if (_dxfEliminateStandaloneVertices(nV,v,new_nV,new_v,nT,*new_t,vlut)) { /* we create a fathers array that is large enough to contain all the triangle corners */ int *fathers = (int *) DXAllocateZero (3 * nT * sizeof (int)); if (fathers) { if (_dxfJoinTriangleCorners(nT, *new_t, e_table, n_edges, edges, fathers)) { if (_dxfBuildOrientedManifold(new_nV, new_v, nT, *new_t, vlut, n_edges, edges, e_table, fathers)) { success = 1;} } DXFree((Pointer)fathers); } } } return success; } /*+--------------------------------------------------------------------------+ | | | _dxfJoinTriangleCorners | | | +--------------------------------------------------------------------------+*/ int _dxfJoinTriangleCorners(int nT, int *t, Htable *e_table, int *n_edges, EdgeS **edges, int *fathers) /* we join two triangle corners only if they share a manifold edge we define a manifold edge as an edge that has at most two incident triangles and such that the two triangles are consistently oriented. Also, the two triangles should not have the same three vertices */ { int success = 0; /* we first determine the edges that are manifold in a second pass, we take all manifold edges and we join the corresponding triangle corners */ register int i; register u_char j,k; int nT3 = 3 * nT; /* allocate an array big enough to contain all possible edges */ EdgeS *edge_array = (EdgeS *) DXAllocateZero (nT3 * sizeof (EdgeS)); if (edge_array) { EdgeS the_edge[1], *cur_edge, *edge_in_table; int *incident_tri = SurfEdgeGetTriangles(the_edge), *endpoints = SurfEdgeGetVertices (the_edge), nE = 0, total_nE = nT3, nE_nonmanifold = 0, /* initialize a pointer to the current triangle */ *the_triangle = t, first_tri, second_tri, corner_v0first, corner_v1first, corner_v0second, corner_v1second, edge_paired_as_manifold = 0, v0, /* surface edge end-points */ v1, *triangles_in_table; /* initialize the fathers array */ _dxfInitFather(nT3, fathers); bzero((char *)e_table, sizeof(Htable)); /* reset the hash table */ /* the size of the table should be chosen such that 2**tbl_size ~ 2 * number of edges for best performances. Here taking nT is plenty enough*/ if (!_dxfInitHashTable(e_table, _dxfFlog2(nT), (int (*)(KEYTYPE))_dxfEdgeSKey, (int (*)(RECTYPE, RECTYPE))_dxfEdgeSFilter, NULL)) goto error; /* initialize the number of faces per edge to 1 */ SurfEdgeLabel(the_edge) = 1; /* for each triangle, add the three edges to the hash table if they are not already indexed in the hash table */ for (the_triangle = t, i=0; i=v1) { endpoints[0] = v0; endpoints[1] = v1; incident_tri[0] = i; incident_tri[1] = -1; } else { endpoints[0] = v1; endpoints[1] = v0; incident_tri[0] = -1; incident_tri[1] = i; } edge_in_table = (EdgeS *)_dxfFindRecord(e_table,(KEYTYPE)the_edge); /* if the edge is not already in the table we add it */ if (!edge_in_table) { /* copy this edge to the edge array */ memcpy(edge_array+nE, the_edge, sizeof(EdgeS)); /* add the edge to the hash table */ if (!_dxfAddRecord(e_table, (RECTYPE)(edge_array+nE), (KEYTYPE)(edge_array+nE))) goto error; /* increment the number of edges */ nE++; } else { /* there is already an edge (v0,v1) in the table */ /* either one of the following two cases apply: CASE ONE: we have already determined that the edge is not a manifold edge, by giving it a label of zero. In that case, we do not even add the edge to the list of edges, but we put a label of zero in the edge slot. and we increment the edge counter by one */ if (SurfEdgeLabel(edge_in_table) == 0) { bzero(edge_array+nE, sizeof(EdgeS)); nE++; } /* CASE TWO: the edge is supposed to be a manifold edge with a label == 2 *but now* there is a third triangle in the picture. so the edge cannot be manifold. so we assign it a label of zero. the problem is, we thought that we had a single edge, and we have two, plus the new one is three. So we increment the edge counter by two. */ else if (SurfEdgeLabel(edge_in_table) == 2) { /* again none of the three edges is a manifold edge */ SurfEdgeLabel(edge_in_table) = 0; bzero(edge_array+nE, sizeof(EdgeS)); nE++; bzero(edge_array+nE, sizeof(EdgeS)); nE++; total_nE++; /*reinstate the edge that we thought was paired */ } /* CASE THREE: we did not establish whether the edge is a manifold edge : label == 1 but we have two triangles and we can establish it */ else if (SurfEdgeLabel(edge_in_table) == 1) { first_tri = i; triangles_in_table = SurfEdgeGetTriangles(edge_in_table); /* if the orientations of the faces do match, we update the edge record to recognize that it is a valid edge, and we will *later* join the corresponding vertices using the father array, if it is established in the edge that the edge (v0,v1) *is indeed* manifold */ if ((v0>=v1) && triangles_in_table[0] == -1) { /* do the triangle orientations match ? */ /* so v0 and v1 are visited in that order in the_triangle. they have to be visited in the opposite order in triangle_in_table */ second_tri = triangles_in_table[1]; /* we know that v1 appears first in second_tri (actually I am not sure what this means because it is only true modulo circular permutation. We can locate v1 in t1: */ j = 0; while (j<3 && t[second_tri*3+j] != v1) j++; /* for the orientation to match, v0 must appear next, right after v1 also, we do not want the two triangles to have exactly the same three vertices */ if (j < 3 && t[second_tri*3+(j+1)%3] == v0 && t[second_tri*3+(j+2)%3] != t[first_tri*3+(k+2)%3]) { edge_paired_as_manifold = 1; triangles_in_table[0] = first_tri; } } else if ((v0=v1) { endpoints[0] = v0; endpoints[1] = v1; incident_tri[0] = i; incident_tri[1] = -1; second_tri = 0; } else { endpoints[0] = v1; endpoints[1] = v0; incident_tri[0] = -1; incident_tri[1] = i; second_tri = 1; } edge_in_table = (EdgeS *)_dxfFindRecord(e_table,(KEYTYPE)the_edge); if (edge_in_table) { triangles_in_table = SurfEdgeGetTriangles(edge_in_table); SurfEdgeLabel(edge_in_table) += 1; triangles_in_table[second_tri] = i; } else { memcpy(edges+nE, the_edge, sizeof(EdgeS)); /* add the edge to the hash table */ if (!_dxfAddRecord(e_table, (RECTYPE)(edges+nE), (KEYTYPE)(edges+nE))) goto error; /* increment the number of edges */ nE++; } } } /* now the final number of edges is known. I return this number to the caller routine */ *n_edges = nE; return 1; error: _dxfFreeHashTable(e_table); if (*new_edges) DXFree((Pointer)(*new_edges)); *new_edges = NULL; return 0; } /*+--------------------------------------------------------------------------+ | | | _dxfEliminateDegenerateTriangles | | | +--------------------------------------------------------------------------+*/ int _dxfEliminateDegenerateTriangles(int nV, int nT, int *t, int *new_nT, int **new_t, int **flut) { int success = 0; /* variable converted to one and returned if the procedure executes successfully */ int trisize = 3 * sizeof (int), *t1 = t, *new_tris = NULL, *face_lut = NULL; register int i; /* initially, set the new number of triangles to zero */ *new_nT = 0; /* 4/20/97: with respect to the original implementation, we now have the difference that there may be a pre-existing "flut" face look-up table, which was created for instance when quads were converted to triangles */ if (*flut == NULL) { face_lut = *flut = (int *) DXAllocateZero(nT * sizeof (int)); if (!(*flut)) goto error; /* initialize the face_lut to point to the identity */ for (i=0;i= 0 && t1[0] < nV && t1[1] >= 0 && t1[1] < nV && t1[2] >= 0 && t1[2] < nV) if (t1[0] != t1[1] && t1[1] != t1[2] && t1[0] != t1[2]) { if (*new_nT < i) face_lut[*new_nT] = face_lut[i]; /* i is always larger (>=) than *new_nT so we will never overwrite a face_lut entry that we will later need */ /* copy this triangle to the new list */ memcpy(new_tris + 3 * (*new_nT), t1, trisize); *new_nT +=1; } /* now reallocate the new triangle array using the new number of triangles reallocate also the face look up table */ if (nT != *new_nT) { *new_t = (int *) DXReAllocate(*new_t, *new_nT * trisize); if (!(*new_t)) goto error; *flut = (int *) DXReAllocate(*flut, *new_nT * sizeof (int)); if (!(*flut)) goto error; } /* don't forget to return "success" if there was no problem */ success = 1; } else goto error; return success; error: if (*new_t) {DXFree((Pointer) (*new_t)); *new_t = NULL;} if (*flut) {DXFree((Pointer) (*flut)); *flut = NULL;} return 0; } /*+--------------------------------------------------------------------------+ | | | _dxfEliminateStandaloneVertices | | | +--------------------------------------------------------------------------+*/ int _dxfEliminateStandaloneVertices(int nV, float *v, int *new_nV, float **new_v_exported, int nT, int *t, int **vlut) { /* the *new_v_exported notation stands for "Manifold set of Vertices" */ register int i, i1; register u_char k; float *new_v = NULL; int *new_tris = t, *vertex_lut = NULL, *old_to_new_vertex_index = NULL; /* allocate the new arrays using the old counts since we potentially eliminate standalone vertices, the arrays cannot grow larger during this procedure (the vertex array may grow larger later on */ int *valence = (int *) DXAllocateZero(nV * sizeof (int)); if (!valence) goto error; old_to_new_vertex_index = (int *) DXAllocateZero(nV * sizeof (int)); if (!old_to_new_vertex_index) goto error; new_v = *new_v_exported = (float *) DXAllocateZero(3 * nV * sizeof (float)); if (!(*new_v_exported)) goto error; vertex_lut = *vlut = (int *) DXAllocateZero(nV *sizeof (int)); if (!(*vlut)) goto error; /* using the set of triangles, compute the valence of each vertex */ _dxfComputeVertexValence(nT, t, valence); *new_nV = 0; /* reset the number of new vertices */ /* build the look-up table of vertices with non-zero valence: the remaining vertices */ for (i=i1=0; i0) { vertex_lut[*new_nV] = i; old_to_new_vertex_index[i] = *new_nV; /* copy the vertices that are not standalone in the new array */ memcpy(new_v + (*new_nV)*3, v + i1, 3 * sizeof (float)); *new_nV += 1; /* increment the count of new vertices */ } /* change the vertex indices in the triangle array */ for (i=i1=0; isize = 1 << numbits; /* the codes will be used as hash entries modulo table->mask */ table->mask = table->size-1; /* alloc the table of entries */ table->bucket = (Nodetype **) DXAllocateZero (sizeof(Nodetype *) * table->size); if (!table->bucket) goto error; for(i=0;isize;i++) table->bucket[i] = (Nodetype *) NULL; table->fcode = fcode; /* function that computes the code from the key */ table->filter = filter;/* filtering function: returns 1 if the filter accepts the record returns 0 if the record is rejected, in which case the next record must be found and processed */ if (compute_key == NULL) /* function that computes the key from the record */ table->compute_key = _dxfIdentifyRecord2Key; else table->compute_key = compute_key; /* initialize the first block of entries */ table->nodes_per_block = 1 << 10 ; table->block = (Block *) NULL; /* initialize the free node */ table->freenode = (Nodetype *) NULL; return 1; error: bzero(table, sizeof (Htable)); return 0; } /*+--------------------------------------------------------------------------+ | | | _dxfNewBlock | | | +--------------------------------------------------------------------------+*/ int _dxfNewBlock(Htable *table) { Block *newblock = NULL; newblock = (Block *) DXAllocateZero (sizeof (Block)); if (!newblock) goto error; newblock->next = table->block; newblock->rec = NULL; table->block = newblock; return 1; error: return 0; } /*+--------------------------------------------------------------------------+ | | | *_dxfgetnode | | | +--------------------------------------------------------------------------+*/ Nodetype *_dxfgetnode(Htable *table) { Nodetype *retval; register u_int i; if (table->freenode) { retval = table->freenode; table->freenode = table->freenode->next; } else { /* create a new block */ /* this can potentially fail if there is not enough memory */ if (!_dxfNewBlock(table)) goto error; /* allocate a new array of nodes */ table->block->rec= (RECTYPE) DXAllocateZero (sizeof(Nodetype) * table->nodes_per_block); if (!table->block->rec) goto error; /* link the nodes newly created */ for(retval= (Nodetype *) table->block->rec, i = 0; i < table->nodes_per_block-1; i++, retval++) { retval->next = retval+1 ; } /* the last node has to point to NULL */ retval->next = (Nodetype *) NULL; retval = (Nodetype *) table->block->rec ; table->freenode = retval+1; } return(retval); error: return (NULL); } /*+--------------------------------------------------------------------------+ | | | *_dxfnewnode | | | +--------------------------------------------------------------------------+*/ Nodetype *_dxfnewnode(Htable *table, RECTYPE rec, int code) { Nodetype *retval = _dxfgetnode(table); /* if _dxfgetnode(table) returns a NULL pointer, we should not try to access it */ if (retval) { retval->code = code; retval->rec = rec; } return(retval); } /*+--------------------------------------------------------------------------+ | | | _dxfHashTableReset | | | +--------------------------------------------------------------------------+*/ void _dxfHashTableReset(Htable *table) { /* resets the hash table by cleaning all buckets and deleting all records */ Block *block=table->block, *next; while(block) { next=block->next; if (block->rec) DXFree((Pointer) block->rec); DXFree((Pointer) block); block=next; } table->block = NULL; /* clear buckets */ if (table->size > 0 && table->bucket) bzero((char *)table->bucket, table->size * sizeof(Nodetype *)) ; return; } /*+--------------------------------------------------------------------------+ | | | _dxfAddRecord | | | +--------------------------------------------------------------------------+*/ int _dxfAddRecord(Htable *table, RECTYPE rec, KEYTYPE key) { Nodetype *node = NULL; register int i; int code= table->fcode(key); node=_dxfnewnode(table,rec,code); if (!node) goto error; /* add the node containing the record to the appropriate bucket */ i=node->code & table->mask; node->next = table->bucket[i]; table->bucket[i]=node; return 1; error: return 0; } /*+--------------------------------------------------------------------------+ | | | *_dxfFindRecord | | | +--------------------------------------------------------------------------+*/ char *_dxfFindRecord(Htable *table, KEYTYPE key) { register int i; table->search_code = table->fcode (key); i = table->search_code & table->mask; table->search_key = key; table->search_node= table->bucket[i]; return (_dxfNextRecord(table)); } /*+--------------------------------------------------------------------------+ | | | *_dxfNextRecord | | | +--------------------------------------------------------------------------+*/ char *_dxfNextRecord(Htable *table) { u_char found=0; Nodetype *node=table->search_node; /* if a filter function was provided, the first node with the same code AND such that the filter function comparing the node's record and the search key "agrees" is selected */ if (table->filter) { while (node && !found) { if (node->code == table->search_code && table->filter(node->rec,table->compute_key(table->search_key))) found=1; else node = node->next; } } /* if no filtering function was provided, the next node with the same code is retrieved */ else while (node && !found) { if (node->code == table->search_code) found=1; else node = node->next; } if (found) { table->search_node = node->next; return ((char *)node->rec); } else return ((char *) NULL); } /*+--------------------------------------------------------------------------+ | | | _dxfIdentifyRecord2Key | | | +--------------------------------------------------------------------------+*/ /* when a filtering function IS provided, we need to compare a current record with the search key. In the general case, a "compute_key" function must be provided to convert the key to a record or something that can be directly compared with a record in many cases the search key is the same as a record, or can be used without conversion. This is when the "dummy" function _dxfIdentifyRecord2Key is used. */ KEYTYPE _dxfIdentifyRecord2Key(RECTYPE rec) { return (KEYTYPE) rec; } /*+--------------------------------------------------------------------------+ | | | _dxfFlog2 | | | +--------------------------------------------------------------------------+*/ /* Log base 2 (integer form) The following relation holds exactly 2 ^ Log2(N) <= N < 2 * 2 ^ Log2(N) if N==0 then the result is by definition -1, but the true result is of course -infty */ int _dxfFlog2(int n) { int d=0; if (n==0) return (-1); while(n>1) {d++; n/=2; } return(d); } /*+--------------------------------------------------------------------------+ | | | _dxfEdgeSKey | | | +--------------------------------------------------------------------------+*/ /* compute the Hash-Table code from the key in the case of a surface edge */ int _dxfEdgeSKey(EdgeS *e) { return(e->v[1]^(e->v[0] << 8)); } /*+--------------------------------------------------------------------------+ | | | _dxfEdgeSFilter | | | +--------------------------------------------------------------------------+*/ int _dxfEdgeSFilter(EdgeS *e1, EdgeS *e2) { if ((e1->v[0] == e2->v[0]) && (e1->v[1] == e2->v[1])) return(1); else return(0); } /*+--------------------------------------------------------------------------+ | | | _dxfEdgeSFilter | | | +--------------------------------------------------------------------------+*/ EdgeS *_dxfFindEdge(Htable *e, int vA, int vB) { EdgeS edge; if (vA > vB) { edge.v[0]=vA;edge.v[1]=vB; } else { edge.v[0]=vB;edge.v[1]=vA; } return (EdgeS *) _dxfFindRecord(e,(KEYTYPE) &edge); } /*+--------------------------------------------------------------------------+ | | | _dxfInitArray | | | +--------------------------------------------------------------------------+*/ void _dxfInitArray(char *array, char *data, int n, int size) { register int i; for(i=0;ivx_data = vertex_data; simp_data->err_volume = old_positional_error; simp_data->edge2index = NULL; simp_data->index2edge = NULL; simp_data->edge_father = NULL; simp_data->normal = NULL; simp_data->area = NULL; simp_data->compactness = NULL; simp_data->valence = NULL; simp_data->boundary_vert = NULL; simp_data->tol_volume = NULL; simp_data->vx_data_error = NULL; simp_data->vx_data_potential_values = NULL; simp_data->edge_heap = NULL; simp_data->edge2index = (int *) DXAllocateZero(simp_data->nE * sizeof (int)); if (!simp_data->edge2index) goto error; simp_data->index2edge = (int *) DXAllocateZero(simp_data->nE * sizeof (int)); if (!simp_data->index2edge) goto error; simp_data->edge_father = (int *) DXAllocateZero(simp_data->nE * sizeof (int)); if (!simp_data->edge_father) goto error; /* initialize the father arrays */ _dxfInitFather(simp_data->nV, simp_data->vertex_father); _dxfInitFather(simp_data->nE, simp_data->edge_father); _dxfInitFather(simp_data->nT, simp_data->triangle_father); simp_data->normal = (Vertex *) DXAllocateZero(simp_data->nT * sizeof (Vertex)); if (!simp_data->normal) goto error; simp_data->area = (float *) DXAllocateZero(simp_data->nT * sizeof (float)); if (!simp_data->area) goto error; simp_data->compactness = (float *) DXAllocateZero(simp_data->nT * sizeof (float)); if (!simp_data->compactness) goto error; /* compute the triangle normals, areas and compactness */ _dxfTrianglesNormalAreaCompactness(simp_data->nT, (Face *)simp_data->triangles, simp_data->vert, simp_data->normal, simp_data->area, simp_data->compactness); /* save the areas in the "old_face_areas" array */ memcpy(simp_data->old_face_areas, simp_data->area, simp_data->nT * sizeof (float)); /* use default values for the maximum angle between face normals and maximum ratio of compactness values */ simp_data->min_scalprod_nor = (float) cos ((double)SIMPLIFY_MAX_ANGLE_NOR); simp_data->compactness_ratio = SIMPLIFY_COMPACTNESS_RATIO; /* allocate and compute the vertex valence */ simp_data->valence = (u_short *) DXAllocateZero(simp_data -> nV * sizeof (u_short)); if (!simp_data->valence) goto error; _dxfComputeVertexValence(simp_data->nT, (int *)simp_data->triangles, simp_data->valence); /* allocate and compute the "boundary_vert" array that flags boundary vertices with a flag value of 1 */ simp_data->boundary_vert = (u_char *) DXAllocateZero(simp_data -> nV * sizeof (u_char)); if (!simp_data->boundary_vert) goto error; _dxfFlagBoundaryVertices(simp_data->nE, simp_data->edge_array, simp_data->boundary_vert); /* allocate and initialize the positional error */ /* there may or may not be positional error defined */ if (simp_data->err_volume) { /* there is positional error already defined */ float *old_error_volume = simp_data->err_volume; simp_data->err_volume = (float *) DXAllocateZero(simp_data->nV * sizeof (float)); if (!simp_data->err_volume) goto error; memcpy(simp_data->err_volume, old_error_volume, simp_data->nV * sizeof (float)); } else { simp_data->err_volume = (float *) DXAllocateZero(simp_data->nV * sizeof (float)); if (!simp_data->err_volume) goto error; } simp_data->tol_volume = (float *) DXAllocateZero(simp_data->nV * sizeof (float)); if (!simp_data->tol_volume) goto error; /* otherwise initialize the tolerance array */ _dxfInitArray((char *)simp_data->tol_volume,(char *) &simp_data->tolerance, simp_data->nV, sizeof (float)); /* treat the vertex-related data */ /* there may or may not be vertex related data defined */ if (simp_data->vx_data && simp_data->data_dim > 0) { register int i; float *old_vertex_data = simp_data->vx_data; simp_data->vx_data = (float *) DXAllocateZero(simp_data->nV * simp_data->data_dim * sizeof (float)); if (!simp_data->vx_data) goto error; simp_data->vx_data_error = (float *) DXAllocateZero(simp_data->nV * sizeof (float)); simp_data->vx_data_potential_values = (float *) DXAllocateZero(3 * simp_data->data_dim * sizeof (float)); if (!simp_data->vx_data_potential_values) goto error; simp_data->vx_old_data = simp_data->vx_data_potential_values + simp_data->data_dim; simp_data->vx_new_data = simp_data->vx_old_data + simp_data->data_dim; /* copy the vertex_data using the vertex_lut */ for (i=0;inV;i++) memcpy(simp_data->vx_data + i * simp_data->data_dim, old_vertex_data + simp_data->vertex_lut[i] * simp_data->data_dim, simp_data->data_dim* sizeof (float)); } else simp_data->vx_data = NULL; /* allocate and initialize the edge heap */ simp_data->edge_heap = _dxfNewHeap(simp_data->nE); if (!simp_data->edge_heap) goto error; simp_data->get_lowest_weight_edge = _dxfRemoveLowestWeightEdgeFromHeap; simp_data->add_edge = _dxfAddEdge2Heap; simp_data->remove_edge = _dxfRemoveEdgeFromHeap; simp_data->heap_initial_index = HeapInitialIndex(simp_data->edge_heap); simp_data->heap_outside_index = HeapOutsideIndex(simp_data->edge_heap); /* initialize the heap indices */ _dxfInitArray((char *)simp_data->edge2index,(char *) &simp_data->heap_initial_index, simp_data->nE, sizeof (int)); /* initialize various counters: */ simp_data->num_edg_collapsed = 0; simp_data->num_edg_weights = 0; simp_data->edg_rejected_4_topology = 0; simp_data->edg_rejected_4_geometry = 0; simp_data->edg_rejected_4_tolerance = 0; /* initialize the last edge added to the heap: */ simp_data->last_edge_added2heap = 0; /* initialize the number of edges left to be tested: */ simp_data->num_edg_remaining_4_testing = simp_data->nE; simp_data->valence_max = SIMPLIFY_VALENCE_MAX; return 1; error: if (simp_data->edge2index) { DXFree((Pointer) simp_data->edge2index); simp_data->edge2index = NULL; } if (simp_data->index2edge) { DXFree((Pointer) simp_data->index2edge); simp_data->index2edge = NULL; } if (simp_data->edge_father) { DXFree((Pointer) simp_data->edge_father); simp_data->edge_father = NULL; } if (simp_data->normal) { DXFree((Pointer) simp_data->normal); simp_data->normal = NULL; } if (simp_data->area) { DXFree((Pointer) simp_data->area); simp_data->area = NULL; } if (simp_data->compactness) { DXFree((Pointer) simp_data->compactness); simp_data->compactness = NULL; } if (simp_data->valence) { DXFree((Pointer) simp_data->valence); simp_data->valence = NULL; } if (simp_data->boundary_vert) { DXFree((Pointer) simp_data->boundary_vert); simp_data->boundary_vert = NULL; } if (simp_data->err_volume && simp_data->err_volume != old_positional_error) { DXFree((Pointer) simp_data->err_volume); simp_data->err_volume = old_positional_error; } if (simp_data->tol_volume) { DXFree((Pointer) simp_data->tol_volume); simp_data->tol_volume = NULL; } if (simp_data->vx_data && simp_data->vx_data != vertex_data) { DXFree((Pointer) simp_data->vx_data); simp_data->vx_data = vertex_data; } if (simp_data->vx_data_error) { DXFree((Pointer) simp_data->vx_data_error); simp_data->vx_data_error = NULL; } if (simp_data->vx_data_potential_values) { DXFree((Pointer)simp_data->vx_data_potential_values); simp_data->vx_data_potential_values = NULL; } if (simp_data->edge_heap) { DXFree((Pointer) simp_data->edge_heap); simp_data->edge_heap = NULL; } return 0; } /*+--------------------------------------------------------------------------+ | | | _dxfFreeSimplificationDataStructure | | | +--------------------------------------------------------------------------+*/ int _dxfFreeSimplificationDataStructure(SimpData *simp_data, float *vertex_data, float *old_positional_error) { if (simp_data->edge2index) { DXFree((Pointer) simp_data->edge2index); simp_data->edge2index = NULL; } if (simp_data->index2edge) { DXFree((Pointer) simp_data->index2edge); simp_data->index2edge = NULL; } if (simp_data->edge_father) { DXFree((Pointer) simp_data->edge_father); simp_data->edge_father = NULL; } if (simp_data->normal) { DXFree((Pointer) simp_data->normal); simp_data->normal = NULL; } if (simp_data->area) { DXFree((Pointer) simp_data->area); simp_data->area = NULL; } if (simp_data->compactness) { DXFree((Pointer) simp_data->compactness); simp_data->compactness = NULL; } if (simp_data->valence) { DXFree((Pointer) simp_data->valence); simp_data->valence = NULL; } if (simp_data->boundary_vert) { DXFree((Pointer) simp_data->boundary_vert); simp_data->boundary_vert = NULL; } if (simp_data->err_volume && simp_data->err_volume != old_positional_error) { DXFree((Pointer) simp_data->err_volume); simp_data->err_volume = old_positional_error; } if (simp_data->tol_volume) { DXFree((Pointer) simp_data->tol_volume); simp_data->tol_volume = NULL; } if (simp_data->vx_data && simp_data->vx_data != vertex_data) { DXFree((Pointer) simp_data->vx_data); simp_data->vx_data = vertex_data; } if (simp_data->vx_data_error) { DXFree((Pointer) simp_data->vx_data_error); simp_data->vx_data_error = NULL; } if (simp_data->vx_data_potential_values) { DXFree((Pointer)simp_data->vx_data_potential_values); simp_data->vx_data_potential_values = NULL; } if (simp_data->edge_heap) { DXFree((Pointer) simp_data->edge_heap); simp_data->edge_heap = NULL; } return 1; } /*+--------------------------------------------------------------------------+ | | | _dxfPartialFreeSimplificationDataStructure | | | +--------------------------------------------------------------------------+*/ int _dxfPartialFreeSimplificationDataStructure(SimpData *simp_data) { if (simp_data->edge2index) { DXFree((Pointer) simp_data->edge2index); simp_data->edge2index = NULL; } if (simp_data->index2edge) { DXFree((Pointer) simp_data->index2edge); simp_data->index2edge = NULL; } if (simp_data->edge_father) { DXFree((Pointer) simp_data->edge_father); simp_data->edge_father = NULL; } if (simp_data->compactness) { DXFree((Pointer) simp_data->compactness); simp_data->compactness = NULL; } if (simp_data->valence) { DXFree((Pointer) simp_data->valence); simp_data->valence = NULL; } if (simp_data->boundary_vert) { DXFree((Pointer) simp_data->boundary_vert); simp_data->boundary_vert = NULL; } if (simp_data->tol_volume) { DXFree((Pointer) simp_data->tol_volume); simp_data->tol_volume = NULL; } if (simp_data->vx_data_error) { DXFree((Pointer) simp_data->vx_data_error); simp_data->vx_data_error = NULL; } if (simp_data->vx_data_potential_values) { DXFree((Pointer)simp_data->vx_data_potential_values); simp_data->vx_data_potential_values = NULL; } if (simp_data->edge_heap) { DXFree((Pointer) simp_data->edge_heap); simp_data->edge_heap = NULL; } return 1; } /*+--------------------------------------------------------------------------+ | | | _dxfTrianglesNormalAreaCompactness | | | +--------------------------------------------------------------------------+*/ int _dxfTrianglesNormalAreaCompactness( int nT, Face *t, Vertex *v, Vertex *t_normal, float *t_area, float *t_compactness) { int success = 1; static VertexD triVert[3]; static double n[4]; register int i; register u_char j,k; /* compute the face normals */ for(i=0;i= 0.) ? 1. : -1., sign2, norm1, norm2; float beta1, beta; static Vertex v1,v2; register u_char k; norm1 = NORM3(x1); memcpy(v1,x1,sizeof(Vertex)); v1[0] += sign1 * norm1; /* we have to record the sign of H1 x1: -sign1 */ if (norm1 >= eps) /* we apply H1 to the vector x2 and replace x2 by H1 x2 if the column is already zero, it is not necessary to use h1 */ { beta1 = - 2. / SCALPROD3(v1,v1); beta = beta1 * SCALPROD3(v1,x2); /* x2 = x2 + beta v1 */ for(k=0;k<3;k++) x2[k] += beta * v1[k]; } /* we create the v2 householder vector; it is by the way only a 2-vector */ {double norm2_d = x2[1]*x2[1] + x2[2] * x2[2]; if (norm2_d > 0.0) norm2 = (float) sqrt(norm2_d); else norm2 = 0.0;} /* we record the sign of H2 x2: -sign2 */ sign2 = (x2[1] >= 0.) ? 1. : -1.; v2[0] = 0; v2[1] = x2[1] + sign2 * norm2; v2[2] = x2[2]; if (norm2 >=eps) { float beta2 = -2. *v2[2] / (v2[1]*v2[1] + v2[2]*v2[2]); n[1] = beta2 * v2[1]; n[2] = 1. + beta2 * v2[2]; } if (norm1 >= eps) { /* we apply H1 to the H2 e3 vector, and that's the normal up to a -1 factor:*/ beta = beta1 * SCALPROD3(v1, n); /* n = n + beta v1 */ for(k=0;k<3;k++) n[k] += beta * v1[k]; } if (sign1 * sign2 <0.) for(k=0;k<3;k++) n[k] = - n[k]; } return 1; } /*+--------------------------------------------------------------------------+ | | | _dxfVectorProductQRD | | | +--------------------------------------------------------------------------+*/ /* compute the vector product using QR decomposition (which, to me, is the same as Householder orthogonalization. I guess the correct name is QR, since a Householder transformation represents a particular symmetry, and several Householder transformation must be combined to produce the final rotation + symmetry. In the case of a 3 * 2 system, either one or two Householder transformations are required. in n[3] store the area of the triangle defined by the origin and the two vectors x1 and x2 */ int _dxfVectorProductQRD(VertexD x1, VertexD x2, double *n) { bzero((char *)n, 4 * sizeof (double)); { static double eps = 1e-15; /* do a QR decomposition of the matrix (x1,x2) */ /* the first householder vector x1 is equal to : x1 + sign(x11) |x1| (1 0 0)^t we will call H1 the corresponding Householder reflexion */ double sign1 = (x1[0] >= 0.) ? 1. : -1., sign2, norm1, norm2; double beta1, beta; static VertexD v1,v2; register u_char k; norm1 = NORM3(x1); memcpy(v1,x1,sizeof(VertexD)); v1[0] += sign1 * norm1; /* we have to record the sign of H1 x1: -sign1 */ if (norm1 >= eps) /* we apply H1 to the vector x2 and replace x2 by H1 x2 if the column is already zero, it is not necessary to use h1 */ { /* beta = -2/ (v1^t v1)* (v1^x2) */ beta1 = - 2. / SCALPROD3(v1,v1); beta = beta1 * SCALPROD3(v1,x2); /* x2 = x2 + beta v1 */ for(k=0;k<3;k++) x2[k] += beta * v1[k]; } /* we create the v2 householder vector; it is by the way only a 2-vector */ norm2 = x2[1]*x2[1] + x2[2] * x2[2]; if (norm2 > 0.0) norm2 = sqrt(norm2); /* we record the sign of H2 x2: -sign2 */ sign2 = (x2[1] >= 0.) ? 1. : -1.; v2[0] = 0; v2[1] = x2[1] + sign2 * norm2; v2[2] = x2[2]; /* compute the area of the triangle */ n[3] = norm1 * norm2 /2.; if (norm2 >=eps) /* we apply H2 to the e3 vector and put the result directly in n */ { /* beta2 = -2 /v2^t v2 * v2^t e3 */ double beta2 = -2. *v2[2] / (v2[1]*v2[1] + v2[2]*v2[2]); n[1] = beta2 * v2[1]; n[2] = 1. + beta2 * v2[2]; } if (norm1 >= eps) { /* we apply H1 to the H2 e3 vector, and that's the normal up to a -1 factor:*/ beta = beta1 * SCALPROD3(v1, n); /* n = n + beta v1 */ for(k=0;k<3;k++) n[k] += beta * v1[k]; } /* There is a unique QR decomposition such that the diagonal of R is made of positive elements and Q is a rotation. The last vector of Q is going to be n. if ( R(1,1)*R(2,2) < 0 ) we have to invert n */ if (sign1 * sign2 <0.) for(k=0;k<3;k++) n[k] = - n[k]; } return 1; } /*+--------------------------------------------------------------------------+ | | | _dxfFlagBoundaryVertices | | | +--------------------------------------------------------------------------+*/ /* as the name indicates, mark the boundary vertices with a value of 1 in the array "boundary_vert" the interior vertices are marked with 0 */ int _dxfFlagBoundaryVertices(int nE, EdgeS *edges, u_char *boundary_vert) { /* loop through edges. if an edge is a boundary edge, then mark both vertices as boundary vertices */ register int i; for (i=0; in && jn) { float key; int pos; key = h->key[i]; h->key[i] = h->key[j]; h->key[j] = key; pos = h->invp[i]; h->invp[i] = h->invp[j]; h->invp[j] = pos; h->perm[h->invp[i]] = i; h->perm[h->invp[j]] = j; } } /*+--------------------------------------------------------------------------+ | | | HeapDown | | | +--------------------------------------------------------------------------+*/ static void HeapDown(int father, Heap *h) { /* trickle an element down the Heap until the Heap condition is verified used when deleting an element from the Heap */ if(0<=father && fatherlast-1) { int left,right,child; while((left=LeftSon(father))last) { child = father; if(h->key[left]key[child]) child = left; if((right=RightSon(father))last && h->key[right]key[child]) child = right; if(child==father) break; HeapSwitch(father,child,h); father = child; } } } /*+--------------------------------------------------------------------------+ | | | HeapUp | | | +--------------------------------------------------------------------------+*/ static void HeapUp(Heap *h) { /* trickle the last heap element up the Heap until the heap condition is verified used when adding an element to the heap */ int child,father; child = h->last-1; while(child>0 && h->key[(father=(child-1)/2)]>=h->key[child]) { HeapSwitch(father,child,h); child = father; } } /*+--------------------------------------------------------------------------+ | | | HeapDel | | | +--------------------------------------------------------------------------+*/ static int HeapDel(int i, Heap *h) { /* delete an element inside a Heap using its actual Heap index this procedure is not meant to be called outside this file */ if (0<=i && ilast) { h->last--; HeapSwitch(i,h->last,h); HeapDown(i,h); return 1; } else return 0; } /*+--------------------------------------------------------------------------+ | | | _dxfHeapDelete | | | +--------------------------------------------------------------------------+*/ int _dxfHeapDelete(int i, Heap *h) { /* delete a Heap element using the original index it was assigned when first added to the Heap. This is the preferred procedure for deleting Heap elements (together with HeapDelMin())*/ if(0<=i && in) return HeapDel(h->perm[i],h); else return 0; } /*+--------------------------------------------------------------------------+ | | | _dxfHeapDelMin | | | +--------------------------------------------------------------------------+*/ int _dxfHeapDelMin(Heap *h) { /* retrieve the Heap element with the smallest key (first element in the Heap), take it off the Heap, and reorganize the Heap */ if(h->last>0) { HeapDel(0,h); return(h->invp[h->last]); } else return(-1); /* heap is empty */ } /*+--------------------------------------------------------------------------+ | | | *_dxfNewHeap | | | +--------------------------------------------------------------------------+*/ Heap *_dxfNewHeap(int n) { /* create a new Heap: allocate and initialize */ Heap *heap = (Heap*)0; if(n>0) { unsigned long bytes = _dxfHeapSize(n); char *block = (char *) DXAllocateZero(bytes); if (!block) goto error; heap = (Heap*)block; heap->n = n; heap->key = (float*)(block+sizeof(Heap)); heap->perm = (int*)(block+sizeof(Heap)+n*sizeof(float)); heap->invp = (int*)(block+sizeof(Heap)+n*sizeof(float)+n*sizeof(int)); _dxfResetHeap(heap); } return(heap); error: return NULL; } /*+--------------------------------------------------------------------------+ | | | _dxfHeapSize | | | +--------------------------------------------------------------------------+*/ unsigned long _dxfHeapSize(int n) { /* used by NewHeap(), compute the size that is needed given the maximum number of elements that can be contained by the Heap */ unsigned long bytes = sizeof(Heap)+n*sizeof(float)+2*n*sizeof(int); return(bytes); } /*+--------------------------------------------------------------------------+ | | | _dxfResetHeap | | | +--------------------------------------------------------------------------+*/ void _dxfResetHeap(Heap *heap) { /* without freeing the memory, re-initialize a Heap for subsequent use */ heap->last = 0; if(heap->n>0) { register int i; for(i=0;in;i++) { heap->key[i]=0.0; heap->perm[i] = heap->invp[i] = i; } } return; } /*+--------------------------------------------------------------------------+ | | | _dxfHeapAdd | | | +--------------------------------------------------------------------------+*/ int _dxfHeapAdd(float key, Heap *h) { int indx = -1; if(h->lastn) /* don't want to enlarge the heap if h->last==h->n ? */ { int last = h->last++; indx = h->invp[last]; h->key[last] = key; HeapUp(h); } return(indx); } /*+--------------------------------------------------------------------------+ | | | _dxfRemoveEdgeFromHeap | | | +--------------------------------------------------------------------------+*/ /* remove a surface edge from the Heap using the edge's original index */ int _dxfRemoveEdgeFromHeap(SimpData *simp_data, int index, int the_edge) { return _dxfHeapDelete(index, simp_data->edge_heap); } /*+--------------------------------------------------------------------------+ | | | _dxfAddEdge2Heap | | | +--------------------------------------------------------------------------+*/ /* add a surface edge to the heap and record the index associated with the edge */ int _dxfAddEdge2Heap(SimpData *simp_data, int the_edge) { int added = 0, nV = simp_data->nV; /* verify that the Heap is not already full and that the edge does not have an index that forbids adding to the Heap */ if (!HeapFull(simp_data->edge_heap) && simp_data->edge2index[the_edge] >=0) { /* compute the representatives of the edge vertices */ int *v = SurfEdgeGetVertices(simp_data->edge_array+the_edge), v0 = FATHER(v[0], simp_data->vertex_father), v1 = FATHER(v[1], simp_data->vertex_father); /* 1) to compute the key to the heap, we choose the edge length to start with */ float weight = DIST3(simp_data->vert[v0], simp_data->vert[v1]); /* 1.1) we add the sum of the errors at the two vertices */ weight += simp_data->err_volume[v0] + simp_data->err_volume[v1]; /* 2) add the edge to the heap */ simp_data->index2edge[ simp_data->edge2index[the_edge]= _dxfHeapAdd( weight, simp_data->edge_heap) ] = the_edge; added = 1; } return added; } /*+--------------------------------------------------------------------------+ | | | _dxfRemoveLowestWeightEdgeFromHeap | | | +--------------------------------------------------------------------------+*/ int _dxfRemoveLowestWeightEdgeFromHeap(SimpData *simp_data) { int idx = _dxfHeapDelMin(simp_data->edge_heap); if (idx != -1) return simp_data->index2edge[idx]; else return -1; } /*+--------------------------------------------------------------------------+ | | | _dxfBuildSimplifiedSurface | | | +--------------------------------------------------------------------------+*/ int _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) /*...........................................................................*/ /* Use the information contained in the vertex and triangle fathers to construct a simplified surface: only the vertices and triangles that are their own "father" are carried to the simplified surface. The vertex_parent and face_parent array are reshuffled to be used later for the purpose of interpolating data (positions-dependent or connections-dependent) from the original surface to the simplified surface) I use indifferently the word "parent" or "father" in this implementation. I think "parent" is better than "father", which refers to the most "direct" parent. But the father of the father is a "parent", and it sounds better to me than a "father". However, there was much existing code using the word "father" so I kept it. */ /*...........................................................................*/ { /* count the new number of vertices */ int *vertex_parents = simp_data->vertex_father, *face_parents = simp_data->triangle_father; register int i,j; float *v = NULL; int *t = NULL, *new_v_number = NULL, *new_t_number = NULL; /* I) build the new vertices and copy the vertex errors and data */ new_v_number = (int *) DXAllocateZero(simp_data->nV * sizeof (int)); if (!new_v_number) goto error; *new_nV = 0; for (i=0; i< simp_data->nV; i++) if (i == FATHER(i,vertex_parents)) { new_v_number[i] = *new_nV; *new_nV += 1; } /* allocate the new vertex array */ v = (float *) DXAllocateZero (*new_nV * sizeof (Vertex)); if (!v) goto error; *new_positional_error = (float *) DXAllocateZero (*new_nV * sizeof (float)); if (!(*new_positional_error)) goto error; if (simp_data->vx_data) { *new_vertex_data = (float *) DXAllocateZero (*new_nV * simp_data->data_dim * sizeof (float)); if (!(*new_vertex_data)) goto error; } /* copy the vertices and other vertex attributes and update the vertex_parents array */ for(j=i=0;inV;i++) if (i == vertex_parents[i]) { memcpy(v + 3 * j,simp_data->vert[i],sizeof(Vertex)); (*new_positional_error)[j] = simp_data->err_volume[i]; if (simp_data->vx_data) memcpy(*new_vertex_data + j * simp_data->data_dim, simp_data->vx_data + i * simp_data->data_dim, simp_data->data_dim * sizeof (float)); j++; } for(i=0;inV;i++) vertex_parents[i] = new_v_number[vertex_parents[i]]; DXFree((Pointer) new_v_number); new_v_number = NULL; /* copy the new vertices into *new_v and reallocate *new_v */ DXFree((Pointer) (*new_v)); *new_v = (float *) DXAllocateZero(*new_nV * sizeof (Vertex)); if (!(*new_v)) goto error; memcpy(*new_v, v, *new_nV * sizeof (Vertex)); DXFree((Pointer) v); v = NULL; /* II) build the new triangles and copy triangle normals and areas */ new_t_number = (int *) DXAllocateZero(simp_data->nT * sizeof (int)); if (!new_t_number) goto error; *new_nT = 0; for (i=0; i< simp_data->nT; i++) if (i == FATHER(i,face_parents)) { new_t_number[i] = *new_nT; *new_nT += 1; } t = (int *) DXAllocateZero (*new_nT * sizeof (Face)); if (!t) goto error; *new_face_areas = (float *) DXAllocateZero (*new_nT * sizeof (float)); if (!(*new_face_areas)) goto error; *face_normals = (float *) DXAllocateZero (*new_nT * sizeof (Vertex)); if (!(*face_normals)) goto error; /* copy the triangles and other triangle attributes and update the face_parents array */ for(j=i=0;inT;i++) if (i == face_parents[i]) { register u_char k; for (k=0;k<3;k++) t[j*3+k] = vertex_parents[ (*new_t)[i * 3 + k]]; (*new_face_areas)[j] = simp_data->area[i]; memcpy(*face_normals + 3 * j, simp_data->normal+i, sizeof(Vertex)); j++; } for(i=0;inT;i++) face_parents[i] = new_t_number[face_parents[i]]; DXFree((Pointer) new_t_number); new_t_number = NULL; /* copy the new triangles into *new_t and reallocate *new_t */ DXFree((Pointer) (*new_t)); *new_t = (int *) DXAllocateZero(*new_nT * sizeof (Face)); if (!(*new_t)) goto error; memcpy(*new_t, t, *new_nT * sizeof (Face)); DXFree((Pointer) t); t = NULL; return 1; error: if (new_t_number) DXFree((Pointer) new_t_number); if (new_v_number) DXFree((Pointer) new_v_number); if (t) DXFree((Pointer) t); if (v) DXFree((Pointer) v); return 0; } /*+--------------------------------------------------------------------------+ | | | _dxfMarkEdgesAdjacent2Boundary | | | +--------------------------------------------------------------------------+*/ /* As the name indicates, I am just using two different flags for edges adjacent to the surface boundary. An edge that is a boundary edge, which is the same as an edge incident to two boundary vertex, is said to be a boundary edge of the second kind, with a flag value of -2 Otherwise, an edge incident to a single boundary vertex is said to be a boundary edge of the first kind, with a flag value of -1. This routine is executed *only* if the user does not wish to simplify the boundary. Once an edge is flagged with a negative flag value, it is impossible to add it later to the edge_heap. */ int _dxfMarkEdgesAdjacent2Boundary(SimpData *simp_data) { int e = 0, nV = simp_data->nV, nE = simp_data->nE, edges_adjacent_2_boundary; u_char *boundary_vertex = simp_data->boundary_vert, boundary1, boundary2; EdgeS *current_edge = simp_data->edge_array; /* for each edge of the surface, determine whether it is a boundary edge of the first or second type. If yes, label it accordingly using the edge2index array */ for (edges_adjacent_2_boundary = nE, e=0; eedge2index[e] = -2; else { int *v = SurfEdgeGetVertices(current_edge); boundary1 = boundary_vertex[v[0]]; boundary2 = boundary_vertex[v[1]]; if (boundary1 || boundary2) simp_data->edge2index[e] = -1; else edges_adjacent_2_boundary--; } } return edges_adjacent_2_boundary; } /*+--------------------------------------------------------------------------+ | | | _dxfBuildEdgeHeap | | | +--------------------------------------------------------------------------+*/ /* add to the Heap edges all the edges that can be added */ int _dxfBuildEdgeHeap(SimpData *simp_data) { int can_add_edge = 1, n_edges_added = 0, nE = simp_data->nE, nV = simp_data->nV; while (can_add_edge && simp_data->last_edge_added2heap < nE) { if ( simp_data->edge2index[simp_data->last_edge_added2heap] == simp_data->heap_initial_index ) /* *new edges* are entered at this stage. if an edge has been visited before, it can only be reentered if a neighbor was simplified. _dxfReinstateEdgesInHeap does this operation */ n_edges_added += (can_add_edge = simp_data->add_edge(simp_data, simp_data->last_edge_added2heap)); simp_data->last_edge_added2heap ++; } /* end while */ return (n_edges_added); } /*+--------------------------------------------------------------------------+ | | | _dxfRotateParentTriangle | | | +--------------------------------------------------------------------------+*/ /* performs a circular permutation (rotation) of the vertex indices of a triangle such that the parent of the first index is the same as "parent_vertex" */ int _dxfRotateParentTriangle(int parent_vertex, int *vertex_fathers, int *tri_v, int *tri_v_parents, int *tri_v_rotated) { /* the triangle must be a parent triangle itself (otherwise all three triangle vertices are not necessarily different */ tri_v_parents[0] = FATHER(tri_v[0], vertex_fathers); tri_v_parents[1] = FATHER(tri_v[1], vertex_fathers); tri_v_parents[2] = FATHER(tri_v[2], vertex_fathers); if (tri_v_parents[0] == parent_vertex) { tri_v_rotated[0] = tri_v[0]; tri_v_rotated[1] = tri_v[1]; tri_v_rotated[2] = tri_v[2]; } else if (tri_v_parents[1] == parent_vertex) { /* rotate the vertex parents */ tri_v_rotated[0] = tri_v_parents[0]; /* use tri_v_rotated[0] as a temporary buffer */ tri_v_parents[0] = tri_v_parents[1]; tri_v_parents[1] = tri_v_parents[2]; tri_v_parents[2] = tri_v_rotated[0]; /* rotate the vertices */ tri_v_rotated[0] = tri_v[1]; tri_v_rotated[1] = tri_v[2]; tri_v_rotated[2] = tri_v[0]; } else { /* rotate the vertex parents */ tri_v_rotated[0] = tri_v_parents[0]; /* use tri_v_rotated[0] as a temporary buffer */ tri_v_parents[0] = tri_v_parents[2]; tri_v_parents[2] = tri_v_parents[1]; tri_v_parents[1] = tri_v_rotated[0]; /* rotate the vertices */ tri_v_rotated[0] = tri_v[2]; tri_v_rotated[1] = tri_v[0]; tri_v_rotated[2] = tri_v[1]; } return 1; } /*+--------------------------------------------------------------------------+ | | | _dxfRotateParentEdge | | | +--------------------------------------------------------------------------+*/ int _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) { /* rotate a parent edge with respect to a parent vertex, so that the child of the parent vertex will be first (index [0]) among the two edge vertices ( currently, by default the first vertex is the vertex with the highest number ) */ edge_v_parents[0] = FATHER(edge_v[0], vertex_fathers); edge_v_parents[1] = FATHER(edge_v[1], vertex_fathers); /* if the vertices are in the right order, just copy the vertices */ if (edge_v_parents[0] == parent_vertex) { edge_v_rotated[0] = edge_v[0]; edge_v_rotated[1] = edge_v[1]; edge_t_rotated[0] = edge_t[0]; edge_t_rotated[1] = edge_t[1]; } else { edge_v_rotated[0] = edge_v_parents[0]; /* use edge_v_rotated[0] as a temporary buffer */ edge_v_parents[0] = edge_v_parents[1]; edge_v_parents[1] = edge_v_rotated[0]; edge_v_rotated[0] = edge_v[1]; edge_v_rotated[1] = edge_v[0]; edge_t_rotated[0] = edge_t[1]; edge_t_rotated[1] = edge_t[0]; } return 1; } /*+--------------------------------------------------------------------------+ | | | _dxfDirectedParentVertexStar | | | +--------------------------------------------------------------------------+*/ /* compute a vertex star and store the triangles in either CLOCKWISE or COUNTER-CLOCKWISE direction */ int _dxfDirectedParentVertexStar(int parent_vertex, int first_triangle, int valence, int *star, SimpData *simp_data, int direction) { /* labeling system for directed triangle stars v*1\ /v*0 e*1\ *1/e*0 *val0-1 \ / *0 ___________\/___________ we omit on purpose the first vertex, because in the case v*val0-1 v0 of a boundary edge, we know it already */ int *vstar = star + valence, *estar = vstar + valence, /* rotated triangle vertices: */ tri_v_rotated[3], tri_v_parents[3], /* rotated edge vertices and adjacent triangles: */ edge_v_rotated[2], edge_v_parents[2], edge_t_rotated[2], star_triangle = first_triangle, nV = simp_data->nV, i; for (i = 0; i < valence; i++) { star[i] = star_triangle; _dxfRotateParentTriangle(parent_vertex, simp_data->vertex_father, simp_data->triangles + 3 * star_triangle, tri_v_parents, tri_v_rotated); /* find previous or next edges and get their parents */ estar[i] = _dxfNextParentEdge(tri_v_rotated, direction); _dxfRotateParentEdge(parent_vertex, simp_data->vertex_father, (int *) SurfEdgeGetVertices (simp_data->edge_array + estar[i]), (int *) SurfEdgeGetTriangles(simp_data->edge_array + estar[i]), edge_v_parents, edge_v_rotated, edge_t_rotated); vstar[i] = edge_v_parents[1]; /* this vertex number should be the same as tri_v_parents[direction] */ if (i < valence-1) /* otherwise, the next triangle would be -1 and finding its father would produce a segmentation fault or another bad thing */ star_triangle = _dxfNextParentTriangle(edge_t_rotated, direction); } return 1; } /*+--------------------------------------------------------------------------+ | | | _dxfManifoldLinkTest | | | +--------------------------------------------------------------------------+*/ /* detect duplicate vertices in the link of an edge. In case duplicate vertices are detected, collapsing an edge would generate a singular vertex. */ int _dxfManifoldLinkTest(int *link0, int *link1, int val0, int val1) { /* among the vertices of the link of v0, link0 and of the link of v1, link1, is there any duplicate vertex ? */ int link_is_manifold = 1, val = val0 + val1; /* since val0 + val1 <= SIMPLIFY_VALENCE_MAX4, I can use a static array here. */ static int link[SIMPLIFY_VALENCE_MAX4]; { register int i; /* gather vertices of link0 and link1 in a single array, sort the array, and determine whether there are any duplicate vertices */ memcpy(link, link0, val0 * sizeof (int)); memcpy(link + val0, link1, val1 * sizeof (int)); qsort(link, val, sizeof(int), (int (*)(const void*, const void*)) _dxfCmpIntSmallFirst); i = 0; while ((i < val -1 ) && (link_is_manifold)) { if (link[i+1] == link[i]) /* if there are duplicate vertices, they must be contiguous in the sorted array */ link_is_manifold = 0; i++; } } return link_is_manifold; } /*+--------------------------------------------------------------------------+ | | | _dxfCmpIntSmallFirst | | | +--------------------------------------------------------------------------+*/ /* compare integers. Using this function in qsort will result in listing the smaller integers first */ int _dxfCmpIntSmallFirst(char *s1, char *s2) { int *f1 = (int *)s1, *f2 = (int *)s2; if (*f2 > *f1) return(-1); else if (*f2 < *f1) return(1); else return(0); } /*+--------------------------------------------------------------------------+ | | | _dxfMakeEdgeBarycenter | | | +--------------------------------------------------------------------------+*/ /* compute the barycenter of two vertices */ void _dxfMakeEdgeBarycenter(Vertex v0, Vertex v1, float alpha0, Vertex v) { register u_char j; float alpha1 = 1. - alpha0; for(j=0; j<3; j++) v[j] = alpha0 * v0[j] + alpha1 * v1[j]; return; } /*+--------------------------------------------------------------------------+ | | | _dxfClosestPointOnEdge | | | +--------------------------------------------------------------------------+*/ /* compute the closest point XP to a point X on an edge (A,B) return the barycentric cordinnates of X with respect to A and B */ float _dxfClosestPointOnEdge(Vertex X, Vertex XP, Vertex A, Vertex B, float *t)/* barycentric coordinates: t, 1-t */ { float d; Vertex BA,BX; float BA_BX,BA_2; MakeVec(B,A,BA); MakeVec(B,X,BX); EdgeBarycentricCoordinates(BA,BX,BA_BX,BA_2); if (BA_BX <= 0.) { *t = 0.; d = NORM3(BX); VecCopy(B,XP); } else if (BA_2 <= BA_BX) { *t = 1.; d = DIST3(A,X); VecCopy(A,XP); } else { /* the following is not robust numerically */ /* d= SCALPROD3(BX,BX) - BA_BX*BA_BX / BA_2; if (d<=0.0) d = 0.0; else d =(float) sqrt ((double) d); */ /* it is more robust to compute d= NORM3(Vecprod(BA,BC))/BA or to just use the distance from X to XP */ *t = BA_BX/BA_2; _dxfMakeEdgeBarycenter(A,B,*t,XP); d = DIST3(X,XP); } return(d); } /*+--------------------------------------------------------------------------+ | | | _dxfNormalize3Vector | | | +--------------------------------------------------------------------------+*/ float _dxfNormalize3Vector(Vertex v) { register u_char j; register double norm; for (norm=0.0,j=0;j<3;j++) norm+=v[j]*v[j]; if (norm>0.0) { norm=sqrt(norm); for (j=0;j<3;j++) v[j] /= norm; } return (float)norm ; } /*+--------------------------------------------------------------------------+ | | | _dxfSolveQuadraticEquation | | | +--------------------------------------------------------------------------+*/ int _dxfSolveQuadraticEquation(double A, double B, double C, double *sol1, double *sol2, double eps) { /* returns the number of real solutions of a second degree equation as well as the solutions everything below eps is supposed to be 0 */ int num_sol; A *= 2.; if ( fabs( A ) < eps) { /* B x + C = 0 */ if ( fabs( B ) < eps) num_sol = 0; else { num_sol = 1; *sol1 = - C / B; } } else { double delta = B * B - 2. * A * C; if ( delta < - eps ) num_sol = 0; else { /* the discriminant is zero or positive */ if (delta < eps) { num_sol = 1; *sol1 = -B / A; } else { delta = sqrt (delta); num_sol = 2; *sol1 = (-B - delta ) /A; *sol2 = (-B + delta ) /A; } } } return num_sol; } /*+--------------------------------------------------------------------------+ | | | _dxfSimplifiedBoundaryVertexLocation | | | +--------------------------------------------------------------------------+*/ /* determine the location of a potential simplified boundary vertex. applies only when a boundary edge is collapsed, which implies that the two endpoints of the edge are boundary vertices */ int _dxfSimplifiedBoundaryVertexLocation(SimpData *simp_data, int v0, int v1, int *vstar0, int *vstar1, int val0, int val1, int method) { int success = 1; if (method == 0) memcpy(simp_data->simplified_vertex, simp_data->vert[v0], sizeof (Vertex)); else if (method == 1) { /* place the simplified vertex so that the length of the boundary will be preserved. The vertex must be situated on an ellipse ___________\/___________\/__________ v*val0-1 v0 v1 v*val1-1 == vl == vr Firstly, we retrieve the vertices vl, v0, v1, vr from the stars. we compute the length of the polygonal arc vl--v0--v1--vr, and call 2a this length The focal points of the ellipsoid are vr and vl, so we note 2c the distance vr-vl. */ float a, a2, c, eps1 = (float)1e-6, *ve = (float *) simp_data->simplified_vertex; Vertex v_zero, v_one, vr, vl, vm; memcpy(v_zero, simp_data->vert[v0], sizeof (Vertex)); memcpy(v_one, simp_data->vert[v1], sizeof (Vertex)); memcpy(vl, simp_data->vert[vstar0[val0-1]], sizeof (Vertex)); memcpy(vr, simp_data->vert[vstar1[val1-1]], sizeof (Vertex)); /* compute the parameters of the ellipsoid */ a = ( DIST3(simp_data->vert[v0], simp_data->vert[v1]) + DIST3(simp_data->vert[vstar0[val0-1]], simp_data->vert[v0]) + DIST3(simp_data->vert[vstar1[val1-1]], simp_data->vert[v1]) ) / 2.; a2 = a * a; c = DIST3(vl,vr)/2.; /* Compute the middle of v_zero and v_one */ vm[0] = (v_zero[0] + v_one[0]) / 2.; vm[1] = (v_zero[1] + v_one[1]) / 2.; vm[2] = (v_zero[2] + v_one[2]) / 2.; /* Then, our first next test is to determine whether a and c are sufficiently different, otherwise, taking the middle of v_zero and v_one will certainly work since we will have the same cord length along the boundary */ if (a == 0.0 || /* if a == 0, the length to preserve is zero, then the segment middle is fine! */ 1. - c/a < eps1) /* limit of the floating point precision */ { memcpy(ve, vm, sizeof(Vertex)); } else /* if (a > c ) */ { float t, x, y, ye; Vertex vp; /* then we consider the middle of v0 and v1, vm, we project vm on the segment vr vl and obtain vp finally, we project vp on the ellipsoid defined by vr, vl, a and c along the direction vp--->vm ve vm |y vl________________0____|_____________vr -c x c vp by projecting vm onto the segment (vl,vc), we obtain vp, and the coordinates of vm in the coordinate system formed by 0(the middle of vr and vl) vr and vm. x and y are the coordinates of vm in this coordinate system We keep x unchanged, and compute ye such that x^2/a^2 + ye^2/b^2 = 1, with b^2 = a^2-c^2 x^2 hence ye^2 = ( 1 - --- ) ( a^2 - c^2 ) a^2 We take for ye the positive square root of ye^2 and finally ve is obtained by ve-vp = ye/y (vm-vp) ve = vp + ye/y (vm-vp) 5/7/97 bug of Floating Point exception noticed on DEC A special case is encountered if c == 0, then we also have x = 0 and the formula is simply ye=a this case is handled with the general case without particular attention, except that we can't divide by c and also we should not execute the part of the code relative to t == 0 or t == 1 */ /* compute x and y of vm with respect to vr and vl */ y = _dxfClosestPointOnEdge(vm, vp, vl, vr, &t); /* what do we do if y == 0, in which case vm is accidentally exactly on the line segment vr-vl? In that case, we still have a non-degenerate ellipse since we ruled out the fact that a and c are the same. Then either v_zero and v_one cannot be on the line segment vr-vl. So we project instead v0 on the segment vr vl v0 ___|________________ vl vp0 vp vr -----> and we use the direction vp0 v0 for projecting the point vp on the ellipse */ if (y <= eps1 * c) /* it can be that c == 0 and y == 0 */ { /* limit of floating point precision */ /* In that case, either v_zero or v_one cannot be on the line vl, vr, so we replace vm with v_zero */ memcpy(vm, v_zero, sizeof(Vertex)); y = _dxfClosestPointOnEdge(vm, vp, vl, vr, &t); if (y <= eps1 *c) { memcpy(vm, v_one, sizeof(Vertex)); y = _dxfClosestPointOnEdge(vm, vp, vl, vr, &t); } } /* what would it mean for y to be equal to zero after these various projections: that all points v_zero, v_one, vm are on vl_vr, meaning that a = c. Just in case, we use an if (y!=0.0) in the end */ /* in case vm does not project orthogonally on the segment (vl,vr) the routine ClosestPointOnEdge projects vm either on vl and assigns 1 to t or on vr and assigns 0 to t we need to determine the angle that vm,vp makes with (vl,vr) */ if (c > 0.0 /* c == 0 is handled as the general case */ && (t == 0. || t == 1.)) { Vertex u,v; double cos_uv, sin_uv, cos_uv_2, sin_uv_2, A, B, C, sol1, sol2, eps2 = 1e-12; MakeVec(vl,vr,u); /* if the segment (vr,vl) is of zero length, which is possible in nasty cases, u will be a vector of zero lenght, and this invalidates the following computations we need to find another solution */ _dxfNormalize3Vector(u); MakeVec(vp,vm,v); _dxfNormalize3Vector(v); cos_uv = u[0] * v[0] + u[1]*v[1] + u[2]*v[2]; /* FScalProd(u,v); */ if (cos_uv < 0.) cos_uv = -cos_uv; cos_uv_2 = cos_uv * cos_uv; sin_uv_2 = 1. - cos_uv_2; /* ye is the solution of the following equation of the second degree: x^2/a^2 + y^2/(a^2-c^2) = 1 <=> (ye cos_uv + c)^2/a^2 + (ye sin_uv)^2/(a^2-c^2) = 1 <=> ye^2( cos_uv_2 + a^2/(a^2-c^2) sin_uv_2) + ye (2 c cos_uv) + c^2-a^2 = 0 A ye^2 + cb ye + C = 0 */ A = cos_uv_2 + sin_uv_2 * a2 / (a+c)/(a-c); B = 2. * cos_uv * c; C = (c+a)*(c-a); /* c^2-a^2 */ _dxfSolveQuadraticEquation(A, B, C, &sol1, &sol2, eps2); /* there is only one positive solution */ if (sol1 > 0.) ye = sol1; else if (sol2 > 0.) ye = sol2; else { ye = 0; /* The quadratic equation has no positive root */ } } /*endif vm project on an extremity */ else /* vm projects normally inside the segment */ { double tmp; /* determine x from t: */ /* vp = t * vl + (1-t) * vr t == 1 <==> x = -c t == 0 <==> x = c x = c - 2 c t = c ( 1-2t) */ x = c * (1. - 2. * t); /* ye = sqrt ((double) (a2 - c * c)*( 1. - x*x/a2)); */ /* since a and c are likely to be very close, a more robust computation is to avoid computations such as a2 - c^2 and compute sqrt(a2) is a is readily available */ tmp = (double) (a -c) * (a+c) * (a-x) * (a+x); /* this is to avoid floating point exceptions when computing the square root */ if (tmp >0.0) ye = sqrt (tmp) / a; else ye = 0.0; } if (y != 0.0) ye /= y; /* do the following vector composition without a function call ve = vp + ye/y (vm-vp) */ ve[0] = vp[0] + ye * (vm[0] - vp[0]); ve[1] = vp[1] + ye * (vm[1] - vp[1]); ve[2] = vp[2] + ye * (vm[2] - vp[2]); /* test whether the boundary length is preserved */ { float l_before = 2. * a, l_after = DIST3(vl,ve) + DIST3(ve,vr), eps3 = simp_data->tolerance * eps1; /* this is a relative accuracy tolerance * floating point accuracy using eps1 here would not be general enough */ if (l_before -l_after > eps3 || l_before - l_after < -eps3) { success = 0; } } } /* if (a>c) otherwise, vm, the middle is good for ve */ } /* if (method = 1, length preservation) otherwise */ return success; } /*+--------------------------------------------------------------------------+ | | | _dxfFastCompactness3DTriangle2 | | | +--------------------------------------------------------------------------+*/ /* compute the triangle compactness and area this is a fast method, in comparison with the method that orthogonalizes the triangle and computes the area afterwards this fast method uses a cross product */ float _dxfFastCompactness3DTriangle2(Vertex *tri, float *triangle_area) { /* compute the squared lenghts of the three edges: */ double sum_lengths_sq = 0., tmp; { register u_char i,j; for(i=0;i<3;i++) { double len_edg_sq = 0; for(j=0;j<3;j++) { tmp = tri[(i+1)%3][j] - tri[i][j]; len_edg_sq += tmp * tmp; } sum_lengths_sq += len_edg_sq; } } { double m0 = tri[1][1] - tri[0][1], m1 = tri[2][2] - tri[1][2], m2 = tri[1][2] - tri[0][2], m3 = tri[2][1] - tri[1][1], m4 = tri[2][0] - tri[1][0], m5 = tri[1][0] - tri[0][0], n0 = m0 * m1 - m2 * m3, n1 = m2 * m4 - m5 * m1, n2 = m5 * m3 - m0 * m4; tmp = n0 * n0 + n1 * n1 + n2 * n2; if (tmp > 0.0) *triangle_area = sqrt (tmp) / 2.; else *triangle_area = 0; } return (float) (FOUR_SQRT_3 * *triangle_area / sum_lengths_sq); } /*+--------------------------------------------------------------------------+ | | | _dxfBoundaryCollapse2KndGeometricallyOk | | | +--------------------------------------------------------------------------+*/ int _dxfBoundaryCollapse2KndGeometricallyOk(SimpData *simp_data, int v0, int v1, /* this procedure must be run twice, once for the star of the vertex v0 and once for the star of the vertex v1 Here, v1 rather stands for the first vertex of the star */ int *star0, int *vstar0, int val0, Vertex *s_normal0, float *s_area0, int direction, float min_scalprod_nor, float compactness_ratio) { int collapsible = 1, i = 0; Vertex faceVert[3]; float *s_comp0 = s_area0 + val0, min_compactness_before, min_compactness_after; /* labeling system for directed triangle stars v*1\ /v*0 e*1\ *1/e*0 *val0-1 \ / *0 ___________\/___________ we omit on purpose the first vertex, because in the case v*val0-1 v0 v1 of a boundary edge, we know it already */ /* retrieve the compactness of the first triangle */ min_compactness_before = simp_data->compactness[star0[0]]; min_compactness_after = 1; /* initialization of the compactness after simplification */ while (collapsible && i < val0-1) { /* retrieve the face compactness before simplification */ min_compactness_before = MIN ( simp_data->compactness[star0[i+1]], min_compactness_before); /* recompute the compactness after changing the first vertex */ memcpy(faceVert[0], simp_data->simplified_vertex, sizeof (Vertex)); if (direction == DIRECTION_CLOCKWISE) { /* reverse the order of triangle vertices, so that the normal that we will compute will respect the surface orientation */ memcpy(faceVert[2], simp_data->vert[vstar0[i]], sizeof (Vertex)); memcpy(faceVert[1], simp_data->vert[vstar0[i+1]], sizeof (Vertex)); } else { memcpy(faceVert[1], simp_data->vert[vstar0[i]], sizeof (Vertex)); memcpy(faceVert[2], simp_data->vert[vstar0[i+1]], sizeof (Vertex)); } s_comp0[i+1] = _dxfFastCompactness3DTriangle2(faceVert, s_area0+i+1); min_compactness_after = MIN ( s_comp0[i+1], min_compactness_after); /* recompute the face normal assuming that the vertex v0 has moved and save this normal for future use*/ _dxfTriangleNormalQR2(faceVert, s_normal0[i+1]); /* check whether the triangle normal orientation is consistent */ collapsible = ((SCALPROD3(s_normal0[i+1], simp_data->normal[star0[i+1]])) > min_scalprod_nor); i++; } if (collapsible) /* check whether the minimum triangle compactness is too much degraded */ collapsible = min_compactness_after > (compactness_ratio * min_compactness_before); return collapsible; } /*+--------------------------------------------------------------------------+ | | | _dxfErrorTestDataOnSurfaceUpdate | | | +--------------------------------------------------------------------------+*/ /* given a new vertex of the mapping between simplified and original surfaces, update the potential data error value at the potential simplified vertex location */ int _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) { int do_collapse = 1; float o_alpha[3], n_alpha[3], *o_data[3], *n_data[3], o_data_error[3], data_difference, new_error_link, tmp_data_diff, new_error_sv, old_error; register int j; register u_char k; o_data[0] = simp_data->vx_data + old_v1 * simp_data->data_dim; o_data[1] = simp_data->vx_data + old_v2 * simp_data->data_dim; o_data[2] = simp_data->vx_data + old_v3 * simp_data->data_dim; n_data[0] = simp_data->vx_data_potential_values; n_data[1] = simp_data->vx_data + new_v2 * simp_data->data_dim; n_data[2] = simp_data->vx_data + new_v3 * simp_data->data_dim; o_alpha[0] = alpha_o_v1; o_alpha[1] = alpha_o_v2; o_alpha[2] = alpha_o_v3; n_alpha[0] = alpha_sv; n_alpha[1] = alpha_n_v2; n_alpha[2] = alpha_n_v3; o_data_error[0] = simp_data->vx_data_error[old_v1]; if (old_v > 1) o_data_error[1] = simp_data->vx_data_error[old_v2]; if (old_v > 2) o_data_error[2] = simp_data->vx_data_error[old_v3]; /* interpolate vertex data in the old and new surfaces and decide if they are sufficiently similar */ for (data_difference = 0., j=0; jdata_dim; j++) { simp_data->vx_old_data[j] = 0.0; for (k=0 ; k < old_v; k++) simp_data->vx_old_data[j] += o_alpha[k] * o_data[k][j]; simp_data->vx_new_data[j] = 0.0; for (k=0 ; k < new_v; k++) simp_data->vx_new_data[j] += n_alpha[k] * n_data[k][j]; tmp_data_diff = simp_data->vx_new_data[j] - simp_data->vx_old_data[j]; data_difference += tmp_data_diff * tmp_data_diff; } /* there are two cases, either the data dimension is one or it is larger than one in the first case, no square root is required for computing the discrepancy between data values, and since I am very stingy about expensive floating point computations, I am avoiding taking this square root if I can. */ if (simp_data->data_dim == 1) {data_difference = (tmp_data_diff > 0.0) ? tmp_data_diff : -tmp_data_diff;} else if (data_difference > 0.0) data_difference = sqrt (data_difference); old_error = 0.0; for (k=0 ; k < old_v; k++) old_error += o_alpha[k] * o_data_error[k]; /* can the sum of the data difference and the old error at the current mesh point be explained by a new error at the simplified vertex, and if it is the case, what is this new error ? */ new_error_link = (new_v == 3) ? alpha_n_v2 * simp_data->vx_data_error[new_v2] + alpha_n_v3 * simp_data->vx_data_error[new_v3]: (new_v == 2) ? alpha_n_v2 * simp_data->vx_data_error[new_v2]: 0.0; /* there is a special case if the simplified vertex is associated with an alpha that's zero: in this case, the error must be explained with the error at the link */ if (alpha_sv == 0. ) { if (new_error_link < old_error + data_difference) do_collapse = 0; /* but otherwise no potential error is computed */ } else /* alpha_sv * error_sv + new_error_link > old_error + data_difference */ { new_error_sv = (old_error + data_difference - new_error_link) / alpha_sv; if (new_error_sv > simp_data->vx_data_potential_err) simp_data->vx_data_potential_err = new_error_sv; do_collapse = (simp_data->vx_data_potential_err <= simp_data->data_tolerance); } return do_collapse; } /*+--------------------------------------------------------------------------+ | | | _dxfSegmentIntersection | | | +--------------------------------------------------------------------------+*/ float _dxfSegmentIntersection(Vertex A, Vertex B, Vertex C, Vertex D, float *lambda, float *mu, int changed_AorB, int changed_C, int changed_D) { /* check whether for two segments in R3, the intersection is in the interior of each of them. lamda values of 1 and 0 qualify as being in the interior, except when one of two edges is of length zero, in which case we set both lambda and mu to be -1, to reflect the fact that there is no intersection */ float dist = 0; static VertexD w[3]; static double work[2], max_z, min_z, sign = 0., scale; register u_char i,j; /* if both Vertices A and B did not change since the previous call of the routine, the scaling factor, householder vector and sign are still correct and do not need to be re-evaluated (this spares us the computation of one square root */ /* if C or D changed, we only have to recompute the portion corresponding to C or D in surface simplification, often the only variable that changes is D */ if (changed_AorB) { for (j=0;j<3;j++) { /* substract A from the rest of the vertices */ w[0][j] = B[j]-A[j]; w[1][j] = C[j]-A[j]; w[2][j] = D[j]-A[j]; } /* record the norm of the first vector and scale the vertices accordingly */ scale = w[0][0] * w[0][0] + w[0][1] * w[0][1] + w[0][2] * w[0][2]; if (scale > 0.0) scale = sqrt (scale); /* particular case if the scale is zero */ if (scale > 0.) { /* apply a scaling factor so that w[0][2] will be negative */ sign = (w[0][2] > 0.) ? -1. : 1.; sign /= scale; for (i=0;i<3;i++) for (j=0;j<3;j++) w[i][j] *= sign; /* apply the scaling */ /* determine the Householder vector for transforming w[0] to z */ w[0][2] -= 1.; _dxfHouseholderPreMultiplicationDbl((double *) w[1], 3, (double *) w[0], 3, 2, work); /* now the first segment is the same as (0,0,1) */ } } else { if (changed_C) { for (j=0;j<3;j++) w[1][j] = (C[j]-A[j]) * sign; /* apply scaling */ _dxfHouseholderPreMultiplicationDbl((double *) w[1], 3, /* apply Householder */ (double *) w[0], 3, 1, work); } if (changed_D) { for (j=0;j<3;j++) w[2][j] = (D[j]-A[j]) * sign; /* apply scaling */ _dxfHouseholderPreMultiplicationDbl((double *) w[2], 3, /* apply Householder */ (double *) w[0], 3, 1, work); } } #define xC_ w[1][0] #define yC_ w[1][1] #define zC_ w[1][2] #define xD_ w[2][0] #define yD_ w[2][1] #define zD_ w[2][2] if (scale > 0.) { /* there are two cases in which there is no intersection if zC and zD >1 or if zC and zD < 0 */ if (zC_> zD_) {max_z = zC_; min_z = zD_;} else {max_z = zD_; min_z = zC_;} if ((min_z > 1.) || (max_z < 0.)) { *lambda = *mu = -1.; } else { static double xDC, yDC, zDC, dlambda, dmu, a, half_b, c, tmp; /* x = xD + mu xDC y = yD + mu yDC */ xDC = xC_ - xD_ ; yDC = yC_ - yD_ ; zDC = zC_ - zD_ ; c = xD_ * xD_ + yD_ * yD_; a = xDC * xDC + yDC * yDC; if (a > 0.) { half_b = xD_ * xDC + yD_ * yDC; dmu = -half_b / a; *mu = dmu; /* mu = -b/(2a) */ tmp = c + half_b * dmu; if (tmp >0.0) tmp = sqrt (tmp); else tmp= 0.0; dist = scale * tmp; /* c - b^2/4a */ /* lambda is given by the 1-z_value at that point */ *lambda = 1. - zD_ - dmu * zDC; } /* otherwise, either the two segment are parallel or the second segment is of zero length */ else if (min_z < max_z) { /* the second segment is parallel to the first *and* there is an intersection */ min_z = MAX(min_z,0.); max_z = MIN(max_z,1.); dlambda = (max_z + min_z)/2.; /* this is definitely between [0,1] and corresponds to a point inside CD to */ *mu = (dlambda - zD_) / zDC; *lambda = 1. - dlambda; dist = scale * sqrt (c); } else { /* the second segment is degenerate */ *lambda = *mu = -1.; } } }/* endif (scale > 0) */ #undef xC_ #undef yC_ #undef zC_ #undef xD_ #undef yD_ #undef zD_ else { /* the first edge is degenerate. There is no mutual intersection */ *lambda = *mu = -1.; } return dist; } /*+--------------------------------------------------------------------------+ | | | _dxfHouseholderPreMultiplication | | | +--------------------------------------------------------------------------+*/ int _dxfHouseholderPreMultiplication( float *A, /* input matrix, given column by column */ int mA, /* leading dimension of matrix a */ float *v, /* householder vector */ int m, /* dimension of v */ int n, /* second dimension of A */ float *w /* workspace, dimension n */ ) { int retval = 0; register int i,j ; float *w1 = w; /* 1) compute beta = -2 vT v: */ float beta = 0; for (i=0;itol_volume[v0], simp_data->tol_volume[v1]); for (i=0;itol_volume[vstar0[i]] < tol_pivot) tol_pivot = simp_data->tol_volume[vstar0[i]]; for (i=1;itol_volume[vstar1[i]] < tol_pivot) tol_pivot = simp_data->tol_volume[vstar1[i]]; /* verify that the maximum tolerance is not already exceeded at the vertices */ i = 0; while (edge_collapsible && i < val0) /* the equality (star_err == tol) is allowed */ if (simp_data->err_volume[vstar0[i++]] > tol_pivot) edge_collapsible = 0; i = 1; while (edge_collapsible && i < val1) if (simp_data->err_volume[vstar1[i++]] > tol_pivot) edge_collapsible = 0; if (simp_data->err_volume[v0] > tol_pivot || simp_data->err_volume[v1] > tol_pivot) edge_collapsible = 0; if (edge_collapsible) { float t, dist; Vertex wp[2]; bzero ((char *)wp, 2 *sizeof (Vertex)); _dxfErrorTestDataOnSurfaceInit(simp_data); /* 1 determine the closest point to the vertex vs on the edge v0,v1 */ dist = _dxfClosestPointOnEdge(simp_data->simplified_vertex, wp[0], simp_data->vert[v0], simp_data->vert[v1],&t); _dxfErrorTestDataOnSurfaceInterpolate2(simp_data, v0, v1, t); potential_err_pivot = dist + t * simp_data->err_volume[v0] + (1.-t) * simp_data->err_volume[v1]; UPDATE_ERR_PIVOT; if (edge_collapsible) { /* 2) determine the closest point to the vertex v0 on the edge (vl,vs) */ dist = _dxfClosestPointOnEdge(simp_data->vert[v0], wp[0], simp_data->vert[vstar0[val0-1]], simp_data->simplified_vertex, &t); /* t * err[vl] + 1-t * err[vs] >= dist + err[v0] */ if (t < 1.) { potential_err_pivot = (dist + simp_data->err_volume[v0] -t * simp_data->err_volume[vstar0[val0-1]] )/ (1.-t); UPDATE_ERR_PIVOT; } /* v0 if (1-t == 0 <==> t = 1) | that means that v0 projects orthogonally on vl vl____________vs we can't treat this constraint and the edge is not collapsible unless err[vl] >= dist + err[v0] */ else if (simp_data->err_volume[vstar0[val0-1]] < dist + simp_data->err_volume[v0]) edge_collapsible = 0; if (edge_collapsible) { if (simp_data->vx_data) edge_collapsible = _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 1, vstar0[val0-1], 0, 1.-t, t, 0., /* vs, vl */ v0, 0, 0, 1., 0., 0.); /* 3) determine the closest point to the vertex v1 on the edge (vs,vr) */ dist = _dxfClosestPointOnEdge(simp_data->vert[v1], wp[0], simp_data->vert[vstar1[val1-1]], simp_data->simplified_vertex, &t); /* t * err[vr] + 1-t * err[vs] >= dist + err[v1] */ if (t < 1.) { potential_err_pivot = (dist + simp_data->err_volume[v1] -t * simp_data->err_volume[vstar1[val1-1]]) / (1.-t); UPDATE_ERR_PIVOT; } else if (simp_data->err_volume[vstar1[val1-1]] < dist + simp_data->err_volume[v1]) edge_collapsible = 0; if (simp_data->vx_data && edge_collapsible) edge_collapsible = _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 1, vstar1[val1-1], 0, 1.-t, t, 0., /* vs, vl */ v1, 0, 0, 1., 0., 0.); if (edge_collapsible) { float lambda1, lambda2; int j; /* 4) take in turn each edge (vs,v*1) (vs,v*val0-1) and intersect them with each edge (v0,v*0), (v0, v*val0-1) */ i = 1; while (edge_collapsible && isimplified_vertex, simp_data->vert[vstar0[i]], simp_data->vert[v0], simp_data->vert[vstar0[j]], &lambda1, &lambda2, !j, 0, 1); /* If the closest point is located in the inside of either edge, determine a crossing point. */ if (lambda2 >= 0. && lambda2 <= 1.) /* inside 2 */ { dist += lambda2 * simp_data->err_volume[v0] + (1. -lambda2) * simp_data->err_volume[vstar0[j]]; if (simp_data->vx_data && lambda1 >= 0. && lambda1 <= 1.) edge_collapsible = _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 2, vstar0[i], 0, lambda1, 1. -lambda1, 0., v0, vstar0[j], 0, lambda2, 1. -lambda2, 0.); if (lambda1 > 0. && lambda1 <= 1.) /* inside 1 */ { potential_err_pivot = (dist - (1. -lambda1) * simp_data->err_volume[vstar0[i]]) / lambda1; UPDATE_ERR_PIVOT; } /* if lambda1 = 0 we can't accomodate that constraint with the error value at the simplified vertex and the edge might be declared "non-collapsible" vs_______v*1 | v0_|_v*0 dist == 0 can occur when the two vertices vstar0[j] and vstar0[i] are the same but we prohibit that */ else if (lambda1 == 0.) { /* the error at vs cannot accomodate that constraint, but maybe it is already satisfied */ if (simp_data->err_volume[vstar0[i]] < dist) edge_collapsible = 0; } } j++; if (i == j) j++; /* we prohibit (i == j) */ } i++; } if (edge_collapsible) { /* 5) do the same as 4) for the star of the vertex v1: take in turn each edge (vs,v*1) (vs,v*val1-1) and intersect them with each edge (v1,v*0), (v1, v*val1-1) */ i = 1; while (edge_collapsible && isimplified_vertex, simp_data->vert[vstar1[i]], simp_data->vert[v1], simp_data->vert[vstar1[j]], &lambda1, &lambda2, !j, 0, 1); if (lambda2 >= 0. && lambda2 <= 1.) /* inside 2 */ { dist += lambda2 * simp_data->err_volume[v1] + (1. -lambda2) * simp_data->err_volume[vstar1[j]]; if (simp_data->vx_data && lambda1 >= 0. && lambda1 <= 1.) edge_collapsible = _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 2, vstar1[i], 0, lambda1, 1. -lambda1, 0., v1, vstar1[j], 0, lambda2, 1. -lambda2, 0.); if (lambda1 > 0. && lambda1 <= 1.) /* inside 1 */ { potential_err_pivot = (dist - (1. -lambda1) * simp_data->err_volume[vstar1[i]]) / lambda1; UPDATE_ERR_PIVOT; } else if (lambda1 == 0.) { /* the error value at vs cannot accomodate that constraint, but maybe it is already satisfied */ if (simp_data->err_volume[vstar1[i]] < dist) edge_collapsible = 0; } } j++; if (i == j) j++; /* we prohibit (i == j) */ } i++; } } } } } } /* update the tolerance and error volumes in case the edge will be simplified */ if (edge_collapsible) { simp_data->err_volume[v0] = largest_err_pivot; /* update the tolerance volume */ simp_data->tol_volume[v0] = MIN(simp_data->tol_volume[v0], simp_data->tol_volume[v1]); _dxfErrorTestDataOnSurfaceEnd(simp_data, v0); } return edge_collapsible; } /*+--------------------------------------------------------------------------+ | | | _dxfRemoveEdgesFromHeap | | | +--------------------------------------------------------------------------+*/ int _dxfRemoveEdgesFromHeap(int *estar, u_short val, int *edge2index, SimpData *simp_data) { int num_removed = 0; register u_short i; /* remove all edges of the star */ for(i=0;iremove_edge(simp_data, idx, the_edge); } return num_removed; } /*+--------------------------------------------------------------------------+ | | | _dxfReinstateEdgesInHeap | | | +--------------------------------------------------------------------------+*/ int _dxfReinstateEdgesInHeap(int *estar, u_short val, SimpData *simp_data, int *num_edg_added) /* return the number of edges that effectively have been reinstated (were in the queue before, were deleted and are now reintered */ { int num_edg_reinstated = 0; register int i; *num_edg_added = 0; /* add all edges of the star to the heap */ for(i=0;iedge2index[estar[i]] == simp_data->heap_outside_index), added = simp_data->add_edge(simp_data, estar[i]); if (added) { *num_edg_added += 1; if (outside) num_edg_reinstated++; } } return num_edg_reinstated; } /*+--------------------------------------------------------------------------+ | | | _dxfCollapseBoundaryEdge2ndKnd | | | +--------------------------------------------------------------------------+*/ int _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) { int n_edges_added[2], n_edges_reinstated = 0; /* e*: v*1\ /v*0 \ / e*1\ *1/e*0 \ *1 / *val0-1 \ / *0 *0 \ / *val1-1 ___________\/___________\/__________ v*val0-1 v0 vs v1 vl vr */ /* father relationships: /\ e*0[0] \e*1[0] / \ /<-----\ / ---\--> *1[1] *1[0] v0 <------- v1 This is only true if val1 > 1 otherwise if val1 ==1 we must have val0 > 1 */ /* 1) update the vertex, triangle and edge fathers */ simp_data->vertex_father[v1] = v0; if (val1 > 1) { simp_data->edge_father [estar1[0]] = simp_data->edge_father [estar0[0]]; simp_data->triangle_father[star1[0]] = simp_data->triangle_father[star1[1]]; } else { simp_data->edge_father [estar0[0]] = simp_data->edge_father [estar1[0]]; simp_data->triangle_father[star0[0]] = simp_data->triangle_father[star0[1]]; } /* 2) change the valences of the vertices */ simp_data->valence[v0] = val0 + val1 - 2; simp_data->valence[vstar0[0]] --; /* 3) remove from the edge heap all edges that have changed */ _dxfRemoveEdgesFromHeap(estar0, val0, simp_data->edge2index, simp_data); _dxfRemoveEdgesFromHeap(estar1, val1, simp_data->edge2index, simp_data); /* mark the 2 discarded edges. verify that estar1[0] was not already counted as outside the heap */ if (simp_data->edge2index[estar1[0]] != simp_data->heap_outside_index) simp_data->num_edg_remaining_4_testing --; simp_data->edge2index[estar1[0]] = -4; simp_data->edge2index[edg_num] = -7; /* 5) reinstate edges in the heap */ n_edges_reinstated = _dxfReinstateEdgesInHeap(estar0, val0, simp_data, n_edges_added); /* reinstate all but the first edge (that collapses) in the star of v1 */ n_edges_reinstated += _dxfReinstateEdgesInHeap(estar1+1, val1-1, simp_data, n_edges_added+1); /* 6) update the number of edges left to be tested */ simp_data->num_edg_remaining_4_testing += n_edges_reinstated; /* 7) update the coordinates of the simplified vertex */ memcpy(simp_data->vert[v0], simp_data->simplified_vertex, sizeof(Vertex)); /* replace the triangle normals, areas, and compactness in the simplified stars */ { register int i; float *s_comp0 = s_area0 + val0, *s_comp1 = s_area1 + val1; /* replace for v0 */ for(i=1;iarea[star0[i]] = s_area0[i]; simp_data->compactness[star0[i]] = s_comp0[i]; memcpy(simp_data->normal[star0[i]], s_normal0[i], sizeof(Vertex)); } /* replace for v1 */ for(i=1;iarea[star1[i]] = s_area1[i]; simp_data->compactness[star1[i]] = s_comp1[i]; memcpy(simp_data->normal[star1[i]], s_normal1[i], sizeof(Vertex)); } } return n_edges_added[0] + n_edges_added[1]; } /*+--------------------------------------------------------------------------+ | | | _dxfCollapsibilityTestsBoundaryEdge2ndKind | | | +--------------------------------------------------------------------------+*/ int _dxfCollapsibilityTestsBoundaryEdge2ndKind(SimpData *simp_data, int edge_num, int v0, int v1, int val0, int val1, int move_simplified_vertex) { int collapsible = 0, *edge_t = (int *) SurfEdgeGetTriangles (simp_data->edge_array+edge_num), *star0 = NULL; Vertex *s_normal0 = NULL; float *s_area0 = NULL; /* perform collapsibility tests particular to a boundary edge of the second kind, such that both edge vertices are boundary vertices In addition of the two vertices touching a boundary, the edge must effectively belong to the boundary, which means that exactly one triangle must be missing, with an index of -1. If both triangles had a -1 index, though, we would not have a valid edge */ if (simp_data->boundary_vert[v0] && simp_data->boundary_vert[v1] && (edge_t[0] == -1 || edge_t[1] == -1)) { /* we only accept the edges such that both end vertices are located on the boundary */ int first_triangle = edge_t[0]; /* perform the valence test. For this particular case, the valence test is 1 <= v(v0) + v(v1) -2 <= valence_max */ if ((val0 + val1 -2 >= 1) && (val0 + val1 -2 < simp_data->valence_max)) { /* build the edge star: contrary to the procedure applied to the surface interior vertices, we compute the vertex star in a directed fashion. Here is the labeling system: v*1\ /v*0 \ / e*1\ *1/e*0 \ *1 / *val0-1 \ / *0 *0 \ / *val1-1 ___________\/___________\/__________ v*val0-1 v0 v1 */ /* allocate storage space for the vertex stars */ star0 = (int *) simp_data->vertex_star_buffer; bzero (star0, (val0 + val1) * 3 * sizeof (int)); /* we next rotate the edge so that edge_t[0] is a valid triangle, t[0] i.e., so that the representation above holds true v[0]-----v[1] t[1] */ if (first_triangle == -1) { /* rotate the edge, i.e. echange v0, v1, val0 and val1 using first_triangle as a temporary storage location assign edge_t[1] to first_triangle */ first_triangle = v0; v0 = v1; v1 = first_triangle; first_triangle = val0; val0 = val1; val1 = first_triangle; first_triangle = edge_t[1]; } {int *vstar0 = star0 + val0, *estar0 = vstar0 + val0, *star1 = estar0 + val0, *vstar1 = star1 + val1, *estar1 = vstar1 + val1; first_triangle = FATHER(first_triangle, simp_data->triangle_father); _dxfDirectedParentVertexStar(v0, first_triangle, val0, star0, simp_data, DIRECTION_COUNTER_CLOCKWISE); /* compute the star of v1 in reverse order */ _dxfDirectedParentVertexStar(v1, first_triangle, val1, star1, simp_data, DIRECTION_CLOCKWISE); /* up to now, we found a parent triangle and we located the next edge, "n_edgeCCW" or "n_edgeCW" depending on whether we are rotating counterclockwise or clockwise around "parent vertex". The parent of this edge must be stored in the e_star. _/ n_edgeCCW/\ parent vertex / \ \_ /____\ n_edge\_ CW Similarly to before, we need to orient the edge so that the parent of the first edge vertex is precisely "parent vertex" */ /* the manifold test is very similar to the test we perform for an interior edge case we compare the vertices of the links and determine whether duplicate vertices exist */ if (_dxfManifoldLinkTest(vstar0+1, vstar1+1, val0-1, val1-1)) { Vertex *s_normal1; float *s_area1 ; s_normal0 = (Vertex *) simp_data->edge_star_buffer_vx; bzero (s_normal0, ( val0 + val1 ) * sizeof (Vertex)); s_normal1 = s_normal0 + val0; s_area0 = (float *) simp_data->edge_star_buffer_fl; bzero (s_area0, 2 * ( val0 + val1 ) * sizeof (float)); s_area1 = s_area0 + 2 * val0; /* position the simplified vertex so that the length of the boundary will be preserved. The vertex must be situated on an ellipse */ if ((_dxfSimplifiedBoundaryVertexLocation(simp_data, v0, v1, vstar0, vstar1, val0, val1, move_simplified_vertex) ) && /* compute the difference in orientation for the triangle normals before and after simplification */ /* compute the minimum compactnesses before and after simplification */ (_dxfBoundaryCollapse2KndGeometricallyOk(simp_data, v0, v1, star0, vstar0, val0, s_normal0, s_area0, DIRECTION_COUNTER_CLOCKWISE, simp_data->min_scalprod_nor, simp_data->compactness_ratio) ) && (_dxfBoundaryCollapse2KndGeometricallyOk(simp_data, v1, v0, star1, vstar1, val1, s_normal1, s_area1, DIRECTION_CLOCKWISE, simp_data->min_scalprod_nor, simp_data->compactness_ratio) ) ) /* geometry test */ { /* compute the tolerance after simplification */ if (_dxfErrorWithinToleranceVBoundary2ndKnd(simp_data, v0, v1, vstar0, vstar1, val0, val1)) { simp_data->num_edg_collapsed +=2; simp_data->num_edg_weights += _dxfCollapseBoundaryEdge2ndKnd(simp_data, edge_num, v0, v1, val0, val1, star0, star1, vstar0, vstar1, estar0, estar1, s_normal0, s_normal1, s_area0, s_area1); collapsible = 1; } else /* reject for tolerance off limits */ simp_data->edg_rejected_4_tolerance ++; } else /* reject for geometry */ simp_data->edg_rejected_4_geometry ++; } else /* reject for topology */ { simp_data->edge2index[edge_num] = -3; /* this condition will not be modified even if the surface is simplified around the edge */ simp_data->edg_rejected_4_topology ++; } } } else simp_data->edg_rejected_4_topology ++; } /* end if (edge_t[0] == -1 || edge_t[1] == -1) */ return collapsible; } /*+--------------------------------------------------------------------------+ | | | _dxfManifoldLinkTest1stKnd | | | +--------------------------------------------------------------------------+*/ int _dxfManifoldLinkTest1stKnd(int v0, int val0, int *star1, int *link1, int val1, SimpData *simp_data) { int link_is_manifold = 1; /* v0 is a boundary vertex v1 is an interior vertex | /\vt |/ \<-----this triangle is *1[0] v0___\_v1 |\ / | \ /<-----this triangle is *1[1] | vb this test works by taking all vertices in the link of v0, except v1, and verify that they are all different from the vertices of the link of v1, except vt, v0 and vb we take from vstar0 as many vertices as val0 (remember that vstar0 has one element more than star0, because v0 is a boundary vertex) we skip the vertices 0, 1, and val1-1 from the link of v1, i.e., we take from vstar1 as many vertices as val1-3 */ int val = val0 + val1 -3; /* since val0 + val1 <= SIMPLIFY_VALENCE_MAX4, I can use a static array here. */ static int link[SIMPLIFY_VALENCE_MAX4]; { int i = 0, tri_v_rotated[3], tri_v_parents[3], /* rotated edge vertices and adjacent triangles: */ edge_v_rotated[2], edge_v_parents[2], edge_t_rotated[2], triangle = star1[0], edge; do { triangle = FATHER (triangle, simp_data->triangle_father); _dxfRotateParentTriangle(v0, simp_data->vertex_father, simp_data->triangles + 3 * triangle, tri_v_parents, tri_v_rotated); edge = _dxfNextParentEdge(tri_v_rotated, DIRECTION_COUNTER_CLOCKWISE); _dxfRotateParentEdge(v0, simp_data->vertex_father, (int *) SurfEdgeGetVertices (simp_data->edge_array + edge), (int *) SurfEdgeGetTriangles(simp_data->edge_array + edge), edge_v_parents, edge_v_rotated, edge_t_rotated); link[i++] = edge_v_parents[1]; triangle = edge_t_rotated[0]; } while (triangle != -1); /* now repeat the same operation in the clockwise direction starting at vb (for Vertex_at_Bottom) */ triangle = star1[1]; do { triangle = FATHER (triangle, simp_data->triangle_father); _dxfRotateParentTriangle(v0, simp_data->vertex_father, simp_data->triangles + 3 * triangle, tri_v_parents, tri_v_rotated); edge = _dxfNextParentEdge(tri_v_rotated, DIRECTION_CLOCKWISE); _dxfRotateParentEdge(v0, simp_data->vertex_father, (int *) SurfEdgeGetVertices (simp_data->edge_array + edge), (int *) SurfEdgeGetTriangles(simp_data->edge_array + edge), edge_v_parents, edge_v_rotated, edge_t_rotated); link[i++] = edge_v_parents[1]; triangle = edge_t_rotated[1]; } while (triangle != -1); /* gather vertices of link0 and link1 in a single array, sort the array, and determine whether there are any duplicate vertices */ memcpy(&link[i], link1+2, (val1 -3) * sizeof (int)); qsort(link, val, sizeof(int), (int (*)(const void*, const void*)) _dxfCmpIntSmallFirst); i = 0; while ((i < val -1 ) && (link_is_manifold)) { if (link[i+1] == link[i]) /* if there are duplicate vertices, they must be contiguous in the sorted array */ link_is_manifold = 0; i++; } } return (link_is_manifold); } /*+--------------------------------------------------------------------------+ | | | _dxfBoundaryCollapse1stKndGeometricallyOk | | | +--------------------------------------------------------------------------+*/ int _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) { int collapsible = 1, i = 1; Vertex faceVert[3]; float *s_comp0 = s_area0 + val0, min_compactness_before, min_compactness_after; /* labeling system *[0]\ / v*[0]____\/____ *[1]/\ / \ *[2] both triangles *[0] and *[1] disappear after edge collapse */ /* compute the compactness of the collapsed triangles */ min_compactness_before = MIN ( simp_data->compactness[star0[1]], simp_data->compactness[star0[0]]); min_compactness_after = 1; /* initialization of the compactness after simplification */ while (collapsible && i < val0-1) { /* retrieve the face compactness before simplification */ min_compactness_before = MIN ( simp_data->compactness[star0[i+1]], min_compactness_before); /* recompute the compactness, and area, after changing the first vertex */ memcpy(faceVert[0], simp_data->simplified_vertex, sizeof (Vertex)); memcpy(faceVert[1], simp_data->vert[vstar0[i]], sizeof (Vertex)); memcpy(faceVert[2], simp_data->vert[vstar0[i+1]], sizeof (Vertex)); s_comp0[i+1] = _dxfFastCompactness3DTriangle2(faceVert, s_area0+i+1); min_compactness_after = MIN ( s_comp0[i+1], min_compactness_after); /* recompute the face normal assuming that the vertex v0 has moved and save this normal for future use*/ _dxfTriangleNormalQR2(faceVert, s_normal0[i+1]); /* check whether the triangle normal orientation is consistent */ collapsible = (SCALPROD3(s_normal0[i+1], simp_data->normal[star0[i+1]]) > min_scalprod_nor); i++; } if (collapsible) /* check whether the minimum triangle compactness is too much degraded */ collapsible = min_compactness_after > (compactness_ratio * min_compactness_before); return collapsible; } /*+--------------------------------------------------------------------------+ | | | _dxfMakeBarycenter | | | +--------------------------------------------------------------------------+*/ void _dxfMakeBarycenter(int nV, Vertex *vv, Vertex bary, float *bary_coord) { register int i; register u_char j; /* initialize barycenter */ for(j=0;j<3;j++) bary[j]=0.0; for(i=0;i 0.) ? 1. : -1, t = sign / ( theta * sign + sqrt (theta * theta + 1.0)), /* we do not check for overflows of theta *theta */ c = 1 / (sqrt (1.0 + t*t)), s = t * c, /* | a b | | c s | | A | | c -s | | | = | | | | | | | b d | | -s c | | C | | s c | */ A = a - t * b, C = d + t * b, Rx = (fabs(A) < eps) ? 0. : (c * e - s * f) / A, Ry = (fabs(C) < eps) ? 0. : (s * e + c * f) / C; *x = c * Rx + s * Ry ; /* apply the inverse rotation */ *y = -s * Rx + c * Ry ; /* " " */ retval = 1; } return retval; } /*+--------------------------------------------------------------------------+ | | | _dxfTriangle3DBarycentricCoordinates2 | | | +--------------------------------------------------------------------------+*/ int _dxfTriangle3DBarycentricCoordinates2(Vertex *tri, Vertex w, Vertex wp, Vertex bary, float *residual) /* compute the barycentric coordinates of a vector with respect to a point the projection wp and the residual are only valid if the function returns 1 (vertex projects inside) */ { int inside; static Vertex x1, x2; static float eps_vertex_outside = EPS_VERTEX_OUTSIDE_TRIANGLE; double a = 0., b = 0., c = 0., d = 0., e = 0., x3, x, y; register u_char j; int origin; /* x1 is the vector representing the shortest edge of the triangle, which is best numerically Min || | x1[0] x2[0] | | w[0] - tri[origin][0] | ||^2 || | x1[1] x2[1] | * | b1 | - | w[1] - tri[origin][1] | || || | x1[2] x2[2] | | b2 | | w[2] - tri[origin][2] | || using the normal equations */ _dxfTriangleBasisVectors(tri, x1, x2, origin); for (j=0; j<3;j++) { x3 = w[j] - tri[origin][j]; a += x1[j] * x1[j]; b += x1[j] * x2[j]; c += x2[j] * x2[j]; d += x1[j] * x3; e += x2[j] * x3; } _dxfSolveSymmetric2x2Eqn(a,b,c,d,e, &x, &y); bary [origin] = 1. - x - y; bary [(origin + 1)%3] = x; bary [(origin + 2)%3] = y; _dxfMakeBarycenter(3, tri, wp, bary); /* determine whether w projects inside or outside tri */ inside = !(bary[0] < eps_vertex_outside || bary[1] < eps_vertex_outside || bary[2] < eps_vertex_outside); /* we actually do not need to create a variable eps_vertex_outside the following test would work all the same inside = !(bary[0] < EPS_VERTEX_OUTSIDE_TRIANGLE || bary[1] < EPS_VERTEX_OUTSIDE_TRIANGLE || bary[2] < EPS_VERTEX_OUTSIDE_TRIANGLE); */ *residual = DIST3(w, wp); return inside; } /*+--------------------------------------------------------------------------+ | | | _dxfClosestPointOnTriangle | | | +--------------------------------------------------------------------------+*/ int _dxfClosestPointOnTriangle(Vertex *tri, Vertex w, Vertex wp, Vertex bary, float *dist) { int inside = _dxfTriangle3DBarycentricCoordinates2(tri, w, wp, bary, dist); if (!inside) { static float eps= EPS_VERTEX_OUTSIDE_TRIANGLE; /* either 2 barycentric coordinates are negative */ if (bary[0] < eps && bary[1] < eps) { bary[0] = bary[1] = 0.; bary[2] = 1.; memcpy(wp, tri[2], sizeof(Vertex)); *dist = DIST3(w,wp); } else if (bary[0] < eps && bary[2] < eps) { bary[0] = bary[2] = 0.; bary[1] = 1.; memcpy(wp, tri[1], sizeof(Vertex)); *dist = DIST3(w,wp); } else if (bary[1] < eps && bary[2] < eps) { bary[1] = bary[2] = 0.; bary[0] = 1.; memcpy(wp, tri[0], sizeof(Vertex)); *dist = DIST3(w,wp); } /* or 1 barycentric coordinate is negative */ else if (bary[0] < eps) { *dist = _dxfClosestPointOnEdge(w, wp, tri[1], tri[2], &bary[1]); bary[2] = 1. - bary[1]; bary[0] = 0.; } else if (bary[1] < eps) { *dist = _dxfClosestPointOnEdge(w, wp, tri[0], tri[2], &bary[0]); bary[2] = 1. - bary[0]; bary[1] = 0.; } else if (bary[2] < eps) { *dist = _dxfClosestPointOnEdge(w, wp, tri[0], tri[1], &bary[0]); bary[1] = 1. - bary[0]; bary[2] = 0.; } } return (inside); } /*+--------------------------------------------------------------------------+ | | | _dxfErrorWithinToleranceVBoundary1stKnd | | | +--------------------------------------------------------------------------+*/ int _dxfErrorWithinToleranceVBoundary1stKnd(SimpData *simp_data, int v0, int v1, int *vstar1, int val1) { int edge_collapsible = 1; register int i; float tol_pivot, potential_err_pivot, largest_err_pivot = 0., dist, lambda1, lambda2; /* set the maximum tolerance for the star */ tol_pivot = simp_data->tol_volume[v1]; for (i=0;itol_volume[vstar1[i]] < tol_pivot) tol_pivot = simp_data->tol_volume[vstar1[i]]; /* verify that the maximum tolerance is not already exceeded at the vertices */ i = 0; while (edge_collapsible && i < val1) if (simp_data->err_volume[vstar1[i++]] > tol_pivot) edge_collapsible = 0; if (simp_data->err_volume[v1] > tol_pivot) edge_collapsible = 0; if (edge_collapsible) { Vertex triangle[3], vs[1],/* local copy of the coordinates of the simplified vertex */ wp[2], alpha[SIMPLIFY_VALENCE_MAX1]; /* val0 + val1 -4 <= valence max <=> val0, val1 <= valence_max +1 because val0 >=3 and val1 >=3 */ float dp; int i_min = -1, link_error_ok = 0; bzero((char *)triangle, 3 * sizeof (Vertex)); bzero((char *)vs, sizeof (Vertex)); bzero((char *)wp, 2 * sizeof (Vertex)); bzero((char *)alpha, SIMPLIFY_VALENCE_MAX1 * sizeof (Vertex)); _dxfErrorTestDataOnSurfaceInit(simp_data); /* the simplified is already positioned, since it must be v0 */ memcpy(vs[0], simp_data->simplified_vertex, sizeof (Vertex)); /* if there is data attached to the surface, we also copy the data of v0 to the new potential data */ if (simp_data->vx_data) memcpy(simp_data->vx_data_potential_values, simp_data->vx_data + v0*simp_data->data_dim, simp_data->data_dim * sizeof(float)); /* 1) determine the closest point to the vertex v1 on the triangles of star1 after simplification */ memcpy(triangle[0], vs[0], sizeof (Vertex)); /* v*[val1-1] labeling system *[0]\ / v*[0]____\/____ *[1]/\ v*[1]/ \ *[2] both triangles *[0] and *[1] disappear after edge collapse */ i = 1; dist = 0.; while (ivert[vstar1[i]], sizeof(Vertex)); memcpy(triangle[2], simp_data->vert[vstar1[i+1]], sizeof(Vertex)); _dxfClosestPointOnTriangle(triangle, simp_data->vert[v1], wp[0], alpha[i], &dp); dp += simp_data->err_volume[v1]; if (alpha[i][0]>0.) { dp = (dp - alpha[i][1] * simp_data->err_volume[vstar1[i]] - alpha[i][2] * simp_data->err_volume[vstar1[i+1]]) / alpha[i][0]; if (i_min == -1 || dp < dist) {dist = dp; i_min = i;} } else { /* otherwise v1 projects onto the link edge of "triangle" the epsilon value at the simplified vertex cannot compensate but it is possible that the error values at the link are already sufficient */ if (alpha[i][1] * simp_data->err_volume[vstar1[i]] + alpha[i][2] * simp_data->err_volume[vstar1[i+1]] >= dp) { /* look no further for the closest triangle */ i_min = i; i = val1-1; /* so that the loop will immediately stop*/ dist = 0; } } i++; } /* the closest triangle wins */ if (i_min >=0) { potential_err_pivot = dist; UPDATE_ERR_PIVOT; if (edge_collapsible && simp_data->vx_data) edge_collapsible = _dxfErrorTestDataOnSurfaceUpdate(simp_data, 3, 1, vstar1[i_min], vstar1[i_min+1], alpha[i_min][0], alpha[i_min][1], alpha[i_min][2], v1, 0, 0, 1., 0., 0.); } else edge_collapsible = 0; /* v1 did not project inside of the new star */ if (edge_collapsible) { /* 2) / v*[val1-1] v*[val1-2] / __/ \ / / __/ \ / v0/__/________ ____\/____ ==\ \__ v*0 /\ vs \ \__ / \ \ \ / \ \ v*1 v*2 take in turn each edge (vs,v*2) ..., (vs, v*[val1-2]) and "intersect" them with each edge (v1, v*[1]) ... (v1, v*[val1-1]) */ int j; i = 2; while (edge_collapsible && i < val1-1) { j = 1; while (edge_collapsible && j < val1) { dist = _dxfSegmentIntersection(simp_data->simplified_vertex, simp_data->vert[vstar1[i]], simp_data->vert[v1], simp_data->vert[vstar1[j]], &lambda1, &lambda2, (j==1), 0, 1); /* If the closest point is located in the inside of either edge, determine a crossing point. */ if (lambda2 >= 0. && lambda2 <= 1.) /* inside 2 */ { dist += lambda2 * simp_data->err_volume[v1] + (1. -lambda2) * simp_data->err_volume[vstar1[j]]; if (simp_data->vx_data && lambda1 >= 0. && lambda1 <= 1.) edge_collapsible = _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 2, vstar1[i], 0, lambda1, 1. -lambda1, 0., v1, vstar1[j], 0, lambda2, 1. -lambda2, 0.); if (lambda1 > 0. && lambda1 <= 1.) /* inside 1 */ { potential_err_pivot = (dist + (lambda1 -1.) * simp_data->err_volume[vstar1[i]]) / lambda1; UPDATE_ERR_PIVOT; } else if (lambda1 == 0.) { /* the error at vs cannot accomodate that constraint, but maybe it is already satisfied */ if (simp_data->err_volume[vstar1[i]] < dist) edge_collapsible = 0; } } j++; if (i == j) j++; /* we prohibit (i == j) */ } i++; } } } if (edge_collapsible) { /* update the tolerance and error values */ /* only the error at the pivot will change */ simp_data->err_volume[v0] = largest_err_pivot; /* update the tolerance volume */ simp_data->tol_volume[v0] = MIN(simp_data->tol_volume[v0], simp_data->tol_volume[v1]); _dxfErrorTestDataOnSurfaceEnd(simp_data, v0); } return edge_collapsible; } /*+--------------------------------------------------------------------------+ | | | _dxfCollapseBoundaryEdge1stKnd | | | +--------------------------------------------------------------------------+*/ int _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) { int n_edges_added = 0, n_edges_reinstated = 0, nV = simp_data->nV, nE = simp_data->nE, et, eb, tri_v_rotated[3], tri_v_parents[3]; /* v*[val1-1] labeling system *[0]\ / v*[0]____\/____ *[1]/\ v*[1]/ \ *[2] both triangles *[0] and *[1] disappear after edge collapse */ /* father relationships: /\ et / \e*1[val1-1] / \ /<-----\ / ---\--> *1[val1-1] *1[0] v0 <------- v1 \ *1[1]/ eb \ --/------> *1[2] \<---/e*1[1] \ / \/ */ /* 1) update the vertex, triangle and edge fathers */ /* 1.1) First retrieve the edges et and eb on "top" and "bottom" of v0, since these edges were not pre-stored */ /* take the triangle star1[0], orient it with respect to v0 and use the oriented triangle vertices to determine the next edge counter clockwise starting from v0: et */ _dxfRotateParentTriangle(v0, simp_data->vertex_father, simp_data->triangles + 3 * star1[0], tri_v_parents, tri_v_rotated); et = (int) (((EdgeS *) _dxfFindEdge(simp_data->e_table, tri_v_rotated[0], tri_v_rotated[2])) - simp_data->edge_array); et = FATHER (et, simp_data->edge_father); /* do the same with the triangle star1[1] to determine the next edge clockwise starting from v0: eb */ _dxfRotateParentTriangle(v0, simp_data->vertex_father, simp_data->triangles + 3 * star1[1], tri_v_parents, tri_v_rotated); eb = (int) (((EdgeS *) _dxfFindEdge(simp_data->e_table, tri_v_rotated[0], tri_v_rotated[1])) - simp_data->edge_array); eb = FATHER (eb, simp_data->edge_father); /* 1.2) */ simp_data->vertex_father[v1] = v0; simp_data->edge_father [estar1[0]] = simp_data->edge_father [estar1[1]] = eb; simp_data->edge_father [estar1[val1-1]] = et; simp_data->triangle_father[star1[0]] = star1[val1-1]; simp_data->triangle_father[star1[1]] = star1[2]; /* 2) change the valences of the vertices */ simp_data->valence[v0] = val0 + val1 - 4; simp_data->valence[vstar1[1]] --; simp_data->valence[vstar1[val1-1]] --; /* 3) remove from the edge heap all edges that have changed */ _dxfRemoveEdgesFromHeap(estar1, val1, simp_data->edge2index, simp_data); /* mark the 2 discarded edges. verify that they were not already counted as outside the heap */ if (simp_data->edge2index[estar1[1]] != simp_data->heap_outside_index) simp_data->num_edg_remaining_4_testing --; if (simp_data->edge2index[estar1[val1-1]] != simp_data->heap_outside_index) simp_data->num_edg_remaining_4_testing --; simp_data->edge2index[estar1[1]] = simp_data->edge2index[estar1[val1-1]] = -4; simp_data->edge2index[edg_num] = -7; /* 5) reinstate all but the two edges and last edge in the star of v1 */ n_edges_reinstated += _dxfReinstateEdgesInHeap(estar1+2, val1-3, simp_data, &n_edges_added); /* 6) update the number of edges left to be tested */ simp_data->num_edg_remaining_4_testing += n_edges_reinstated; /* 7) update the coordinates of the simplified vertex" no need since its coordinates do not change*/ /* replace the triangle normals, areas and compactness in the simplified star of v1 */ { register int i; float *s_comp1 = s_area1 + val1; /* replace for v1 */ for(i=2;iarea[star1[i]] = s_area1[i]; simp_data->compactness[star1[i]] = s_comp1[i]; memcpy(simp_data->normal[star1[i]], s_normal1[i], sizeof(Vertex)); } } return n_edges_added; } /*+--------------------------------------------------------------------------+ | | | _dxfCollapsibilityTestsBoundaryEdge1stKind | | | +--------------------------------------------------------------------------+*/ int _dxfCollapsibilityTestsBoundaryEdge1stKind(SimpData *simp_data, int edge_num, int v0, int v1, int val0, int val1) { /* perform collapsibility tests particular to a boundary edge of the first kind, such that exactly one of the edge vertices is a boundary vertex */ int collapsible = 0; int *star1 = NULL; Vertex *s_normal1 = NULL; float *s_area1 = NULL; /* the valence of v0 becomes val0 + val1 - 4 */ if ((val0 + val1 -4 >= 1) && (val0 + val1 -4 < simp_data->valence_max)) { int *edge_t = (int *) SurfEdgeGetTriangles (simp_data->edge_array+edge_num), nT = simp_data->nT, first_tri = edge_t[0], second_tri = edge_t[1], /* this is the counter-clockwise order with respect to v1 */ *vstar1, *estar1; /* permute v0 and v1 so that v0 is the boundary vertex */ if (simp_data->boundary_vert[v1]) { int tmp = v0; v0 = v1; v1 = tmp; tmp = val0; val0 = val1; val1 = tmp; first_tri = edge_t[1]; second_tri = edge_t[0]; } /* initialize the storage pointer for the buffer star1 */ star1 = (int *) simp_data->vertex_star_buffer; bzero (star1, 3 * val1 * sizeof (int)); /* | / first_triangle |/ v0----------------v1 \ \ boundary */ estar1 = star1 + val1; vstar1 = estar1 + val1; first_tri = FATHER(first_tri, simp_data->triangle_father); /* first star elements (index 0)*/ star1[0] = first_tri; estar1[0] = edge_num; vstar1[0] = v0; /* complete the star from indices 1 to val1-1 */ _dxfParentVertexStar(v1, second_tri, (u_short) val1, star1, simp_data); if (_dxfManifoldLinkTest1stKnd(v0, val0, star1, vstar1, val1, simp_data)) { /* the simplified vertex is already positioned, since it must be v0 */ memcpy(simp_data->simplified_vertex, simp_data->vert[v0], sizeof (Vertex)); /* initialize the storage pointer for the buffer s_normal1 */ s_normal1 = (Vertex *) simp_data->edge_star_buffer_vx; bzero(s_normal1, val1 * sizeof (Vertex)); /* initialize the storage pointer for the buffer s_area1 */ s_area1 = (float *) simp_data->edge_star_buffer_fl; bzero (s_area1, 2 * val1 * sizeof (float)); if (_dxfBoundaryCollapse1stKndGeometricallyOk(simp_data, v1, star1, vstar1, val1, s_normal1, s_area1, simp_data->min_scalprod_nor, simp_data->compactness_ratio)) { /* perform the error tolerance test */ if (_dxfErrorWithinToleranceVBoundary1stKnd(simp_data, v0, v1, vstar1, val1)) { simp_data->num_edg_collapsed +=3; simp_data->num_edg_weights += _dxfCollapseBoundaryEdge1stKnd(simp_data, edge_num, v0, v1, val0, val1, star1, vstar1, estar1, s_normal1, s_area1); collapsible = 1; } else /* reject for tolerance off limits */ simp_data->edg_rejected_4_tolerance ++; } else /* reject for geometry */ simp_data->edg_rejected_4_geometry ++; } else { /* this condition will not change if the surface is further simplified */ simp_data->edge2index[edge_num] = -3; simp_data->edg_rejected_4_topology ++; } } else simp_data->edg_rejected_4_topology ++; return collapsible; } /*+--------------------------------------------------------------------------+ | | | _dxfCollapseTopologicallyFeasible | | | +--------------------------------------------------------------------------+*/ int _dxfCollapseTopologicallyFeasible(SimpData *simp_data, int *vstar0, int *vstar1, u_short val0, u_short val1) { register u_short i, j; int feasible = 1; /* count each vertex only once *1[val1-1] /\ *0[1] counted twice *1[0] / \ *0[0] \ / *1[1] \/ *0[val0-1] counted twice */ { /* verify that no two vertices of the link are the same */ val0 -= 2; val1 -= 2; i = 2; while(( i <= val0) && (feasible)) { j = 2; while ((j <= val1) && (feasible)) { if (vstar0[i] == vstar1[j]) feasible = 0; j++; } i++; } } return(feasible); } /*+--------------------------------------------------------------------------+ | | | _dxfBuildParentEdgeStars | | | +--------------------------------------------------------------------------+*/ void _dxfBuildParentEdgeStars(EdgeS *edg, int v0, int v1, u_short val0, u_short val1, int *star0, int *star1, SimpData *simp_data) /*...........................................................................*/ /* The following procedure computes the star of v0 and v1 belonging to an "active" surface data structure, i.e. the surface is undergoing a simplification process. the ordering of the stars is done such that: e2 e1 e5 e4 \ t2 / t1 \ t5 / t3 \ / t0 \ / t4 e3 _______\/_______e0 e0 _______\/_______e3 /\v0 v1 v0 /\v1 t4 / \ t0 t1 / \ t3 / t5 \ / t2 \ e4 e5 e1 e2 */ /*...........................................................................*/ { int *estar0 = star0 + val0, *vstar0 = estar0 + val0; int *estar1 = star1 + val1, *vstar1 = estar1 + val1, e_num = (int) (edg - simp_data->edge_array); int *edg_triangles = SurfEdgeGetTriangles(edg); /* complete first the straighforward portion of the stars /\ / \ / \ father(edg_triangles[0]); / \ v0/edg->num\ v1 \ / \ / \ / father(edg_triangles[1]); \ / \/ */ star0 [0] = FATHER(edg_triangles[1], simp_data->triangle_father); vstar0[0] = v1; estar0[0] = e_num; /* the simplifiable edge is always a root edge */ _dxfParentVertexStar(v0, edg_triangles[0], val0, star0, simp_data); star1 [0] = FATHER(edg_triangles[0], simp_data->triangle_father); vstar1[0] = v0; estar1[0] = e_num; _dxfParentVertexStar(v1, edg_triangles[1], val1, star1, simp_data); return; } /*+--------------------------------------------------------------------------+ | | | _dxfParentVertexStar | | | +--------------------------------------------------------------------------+*/ void _dxfParentVertexStar(int vf, int t0, u_short val, int *star, SimpData *simp_data) /*...........................................................................*/ /* an "active" vertex is a vertex belonging to the inside of a surface ( not a boundary vertex ) undergoing a simplification process */ /*...........................................................................*/ { int *estar = star + val; int *vstar = estar + val; register int i; int t = t0; /* complete the star from 1 to the valence, the 0-elements have been already completed */ for(i=1;itriangle_father); /* complete the triangle part of the vertex star */ star[i] = t; do { vA = simp_data->triangles[3 * t + j]; vB = simp_data->triangles[3 * t + (j+2)%3]; vAf = FATHER(vA, simp_data->vertex_father); } while((j++ < 3) && (vAf != vf)); /* now complete the edge portion of the vertex star */ edg = _dxfFindEdge(simp_data->e_table,vA,vB); e_num = (int) (edg - simp_data->edge_array); /* differential pointer is also the edge number */ ef = FATHER (e_num, simp_data->edge_father); estar[i] = ef; /* find out the parents of the vertices of ef */ edgf = simp_data->edge_array+ef; endpoints = SurfEdgeGetVertices(edgf); facing_tris = SurfEdgeGetTriangles(edgf); vAf = FATHER(endpoints[0], simp_data->vertex_father); vBf = FATHER(endpoints[1], simp_data->vertex_father); /* and the vertex portion of the vertex star */ if (vAf == vf) { vstar[i] = vBf; /* /\f0 ^ccw direction the next triangle is v0/__\v1 | v0 = vf, v0 > v1 \ / | \/f1 */ t = facing_tris[0]; } else { vstar[i] = vAf; t = facing_tris[1]; } } return; } static void CentroidEdgeLink(Vertex centroid, Vertex *vert, int val0, int val02, int val014); /*+--------------------------------------------------------------------------+ | | | CentroidEdgeLink | | | +--------------------------------------------------------------------------+*/ static void CentroidEdgeLink(Vertex centroid, Vertex *vert, int val0, int val02, int val014) { register int i; register u_char j; bzero((char *)centroid, sizeof(Vertex)); for(j=0;j<3;j++) { for(i=0; i<=val02; i++) centroid[j] += vert[i][j]; for(i=val0; i<=val014; i++) centroid[j] += vert[i][j]; centroid[j] /= val014; } return; } static void BuildEdgeStarVal1Is3(Vertex *, Face *, int, int, int, int, int, int *, Vertex *); /*+--------------------------------------------------------------------------+ | | | BuildEdgeStarVal1Is3 | | | +--------------------------------------------------------------------------+*/ static void BuildEdgeStarVal1Is3(Vertex *vert, Face *face, int v0, int v1, int val0, int val01, int val02, int *vstar0, Vertex *surfvert) { register int i,i1; /* example of numbering for the case when val0 = 6, val1 = 3 v1 _v0 \ _//| \ t1 _/ / | The interior vertices come last t2 \ _/t0/ | The indices of triangles are shifted by one v2______\/___/ t6| with respect to the vertex star v5/\_ \v6 | t3 / \t5\ | / t4 \_\ | / \\| v3 v4 */ /* 1) copy the locations of the vertices */ /* 1.1) star of v0 */ for (i=0, i1=1; i< val01; i++, i1++) memcpy(vert[i], surfvert[vstar0[i1]], sizeof(Vertex)); memcpy(vert[val01], surfvert[v0], sizeof(Vertex)); /* 1.2) star of v1 */ memcpy(vert[val0], surfvert[v1], sizeof(Vertex)); /* 2) create the faces or connections between vertices */ face[0][0] = val01; face[0][1] = val0; face[0][2] = 0; for (i = 1, i1 = 0; i <= val02; i++, i1++) { face[i][0] = val01; face[i][1] = i1; face[i][2] = i; } face[val01][0] = val01; face[val01][1] = val02; face[val01][2] = val0; face[val0][0] = val0; face[val0][1] = val02; face[val0][2] = 0; return; } /*+--------------------------------------------------------------------------+ | | | _dxfBuildEdgeStar | | | +--------------------------------------------------------------------------+*/ void _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) /*...........................................................................*/ /* represent the edge star with origin at the centroid of the edge link copy the vertices and faces of both vertex stars in a single array: here is the new numbering of the vertices: 1 0 val0+val1-4 \ / \ / .. \ / \ / .. ____\/____ ____\/____val0+1 val0-1/\v0 val0+val1-3/\ / \ / \ / \ / \val0 .. val0 -2 The number of triangles in star_val is val0 + val1 -2 The number of vertices is val0 + val1 -4 + 2, since v0 and v1 belong to the edge star In a vertex star, there is one more vertex (the center) than there are triangles. */ /*...........................................................................*/ { /* there are two cases: either val1 is exactly 3 or it is > 3 */ int star_val = val0 + val1 - 2, val013 = star_val -1, val014 = val013 -1, val01 = val0 -1, val02 = val01 -1, *vstar0 = star0 + 2 * val0, *vstar1 = star1 + 2 * val1; Vertex translation; if (val1 == 3) { BuildEdgeStarVal1Is3(star_vert, star_face, v0, v1, val0, val01, val02, vstar0, s->vert); } else { register int i,i1; /* 0) copy the locations of the vertices */ /* 0.1) star of v0 */ for (i=0, i1=1; i< val01; i++, i1++) memcpy(star_vert[i], s->vert[vstar0[i1]], sizeof(Vertex)); memcpy(star_vert[val01], s->vert[v0], sizeof(Vertex)); /* 0.2) star of v1 */ for (i=val0, i1=2; i<= val014; i++, i1++) memcpy(star_vert[i], s->vert[vstar1[i1]], sizeof(Vertex)); memcpy(star_vert[val013], s->vert[v1], sizeof(Vertex)); /* 1) copy the faces */ /* numbering of the faces : val0 + \ 1 / \val1/ 2 \ / 0 \-3/ ____\/____ ____\/____ /\ /\ .. / \ val0-1 / \val0+1 / \ /val0\ */ star_face[0][0] = val01; star_face[0][1] = val013; star_face[0][2] = 0; for (i = 1, i1 = 0; i <= val02; i++, i1++) { star_face[i][0] = val01; star_face[i][1] = i1; star_face[i][2] = i; } star_face[val01][0] = val01; star_face[val01][1] = val02; star_face[val01][2] = val013; star_face[val0][0] = val013; star_face[val0][1] = val02; star_face[val0][2] = val0; for (i = val0+1, i1 = val0; i <= val014; i++, i1++) { star_face[i][0] = val013; star_face[i][1] = i1; star_face[i][2] = i; } star_face[val013][0] = val013; star_face[val013][1] = val014; star_face[val013][2] = 0; } /* 2) copy the normals of the faces and plug them into the equations of the planes */ _dxfCopyFaceNormalsAreasCompactness(star_plane, star_areas, star_comp, star0, star1, val0, val013, s->normal, s->area, s->compactness); /* 3) compute the centroid of the edge link, and put it in the end of the star_vert array (the vertices v0 and v1 are excluded) */ CentroidEdgeLink(star_vert[star_val], star_vert, val0, val02, val014); /* 4) change the coordinate system such that the origin will be the centroid of the edge link */ _dxfOppositeVector(star_vert[star_val], translation); _dxfApplyTranslation(star_val, star_vert, translation); return; } /*+--------------------------------------------------------------------------+ | | | _dxfCopyFaceNormalsAreasCompactness | | | +--------------------------------------------------------------------------+*/ void _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) { register int i, i1; int val01 = val0-1; /* copy the normals of the faces and plug them into the equations of the planes */ for (i=0, i1=1; i< val01; i++, i1++) { memcpy(plane[i], t_normal[star0[i1]], sizeof(Vertex)); s_area[i] = t_area[star0[i1]]; s_comp[i] = t_comp[star0[i1]]; } memcpy(plane[val01], t_normal[star0[0]], sizeof(Vertex)); s_area[val01] = t_area[star0[0]]; s_comp[val01] = t_comp[star0[0]]; for (i=val0, i1=2; i<= val013; i++, i1++) { memcpy(plane[i], t_normal[star1[i1]], sizeof(Vertex)); s_area[i] = t_area[star1[i1]]; s_comp[i] = t_comp[star1[i1]]; } return; } /*+--------------------------------------------------------------------------+ | | | _dxfApplyTranslation | | | +--------------------------------------------------------------------------+*/ void _dxfApplyTranslation(int nV, Vertex *vv, Vertex T) { register int i; register u_char j; for(i=0;i 0.0) for (i=0; i < 4; i++) p[i] /= norm; } return; } /*+--------------------------------------------------------------------------+ | | | _dxfSimplifiedVertexLocation | | | +--------------------------------------------------------------------------+*/ void _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, /* workspace needed for the householder transform: */ float *work, u_char method) /* various methods for computing the position of the simplified vertex */ { if (method == 0) /* v0 is the simplified vertex */ memcpy(s->simplified_vertex, star_vert[val0-1], sizeof(Vertex)); else if (method == 1) { /* take the vertex that minimizes the sum of distances to the vertices of the edge star plus the planes defined by the triangles of the edge star and that preserves the volume of the edge star */ float star_volume, /* volume of the edge star (x 6) */ sign; int star_val_1 = star_val + 1, size_planes = (star_val_1) * sizeof (Plane); float *star_plane_sav = (float *) (star_plane + star_val_1); static Vertex v, z = {0,0,1.}; star_volume = _dxfStarPlaneEquationsAndVolume(star_vert, star_face, star_plane, star_areas, star_val); /* compute the equation of the plane p such that the star of any vertex on that plane would have the same volume as star_volume */ _dxfFPlaneForIdenticalVolume(val0, val1, star_plane[star_val], star_vert, star_face, star_volume); /* now change again the coordinate system such that either (0,0,1) or (0,0,-1) corresponds to the normal of P, change the equations of the planes. THE VERTICES ARE LEFT UNCHANGED */ /* compute the householder vector */ sign = (star_plane[star_val][2] > 0.) ? 1. : -1.; _dxfComposeVectors(star_plane[star_val], z, sign, v); /* Apply the Householder transform to the planes while saving their previous positions */ memcpy(star_plane_sav, star_plane, size_planes); _dxfHouseholderPreMultiplication((float *)star_plane, 4, v, 3, star_val_1, work); /* this version of the positioning solves the Normal Equations rather than the LS equations and has a closed form solution (does not require usage of a numerical library) */ _dxfPositionVertexLSConstraints2( star_plane, star_val, s->simplified_vertex); /* change back the position of the simplified vertex to the original coordinate system */ _dxfHouseholderPreMultiplication(s->simplified_vertex, 3, v, 3, 1, work); /* idem with the planes */ memcpy(star_plane, star_plane_sav, size_planes); } else if (method == 2) { /* take the middle of v0 and v1 for the position of the simplified vertex */ u_short val01 = val0-1, val013 = val0+val1-3; register u_char j; for (j=0;j<3;j++) s->simplified_vertex[j] = (star_vert[val01 ][j] + star_vert[val013][j] ) /2.; } return; } /*+--------------------------------------------------------------------------+ | | | _dxfCollapseGeometricallyFeasible | | | +--------------------------------------------------------------------------+*/ int _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) { int feasible = 1; /* check whether any triangle normal might have switched its orientation */ if (_dxfNoFlipOverCheck(1, val0-1, snor0, splane, min_scalprod) == 1 && _dxfNoFlipOverCheck(val0, sval, snor1, splane, min_scalprod) == 1) { float old_comp_min = _dxfMinFloatArray(old_comp,sval); feasible = (new_comp_min > (compactness_ratio * old_comp_min)); } else feasible = 0; return feasible; } /*+--------------------------------------------------------------------------+ | | | _dxfNoFlipOverCheck | | | +--------------------------------------------------------------------------+*/ int _dxfNoFlipOverCheck(u_short first_face, u_short last_face, Vertex *nor, Plane *plane, float min_scalprod_nor) { register u_short i1=0, i=first_face; int orientation_consistent; do { orientation_consistent = (SCALPROD3(nor[i1],plane[i]) > min_scalprod_nor); i++; i1++; } while ((orientation_consistent ) && (i < last_face)); return orientation_consistent; } /*+--------------------------------------------------------------------------+ | | | _dxfErrorWithinToleranceV | | | +--------------------------------------------------------------------------+*/ int _dxfErrorWithinToleranceV(SimpData *simp_data, Vertex *star_vert, Face *star_face, Plane *star_plane, int *star0, int *star1, int val0, int val1) { /* The purpose of Method V is to construct in parallel a tiling of e* (edge undergoing collapsibility tests) denoted til(e*) and a tiling of vs* (potential simplified vertex) denoted til(vs*), such that there is a one to one mapping between til(e*) and til(vs*) */ /* refresher: numbering strategies for e*: numbering of the vertices : ~~~~~~~~~~~~~~~~~~~~~~~~~ 1 0 val0+val1-4 \ / \ / .. \ / \ / .. ____\/____ ____\/____val0+1 val0-1 /\ val0+val1-3/\ / \ / \ / \ / \ .. val0-2 val0 numbering of the faces : t0,t1,...t_val0+val1-3 triangles of e* ~~~~~~~~~~~~~~~~~~~~~~ val0 + \ 1 / \val1/ 2 \ / 0 \-3/ ____\/____ ____\/____ /\ /\ .. / \ val0-1 / \val0+1 / \ /val0\ q1,...,qval0-2,....,qval0,...,q_val0+val1-3 triangles of vs* */ /* 1) determine the closest point to the vertex vs on the triangles t0 and t_val0-1 */ /* 2) determine the closest point to the vertex v0 on the triangles q1,...,qval0-2. */ /* 3) determine the closest point to the vertex v1 on the triangles qval0,...,q_val0+val1-3 */ /* 4) Take in turn each edge (vs,v0),(vs,v1),..., (vs,v_val0-2) call the ith edge ei determine the closest point on each edge incident to v_val0-1: (v_val0-1, v0),..., (v_val0-1, v_val0-2), and (v_val0-1, v_val0+val1-3) If the closest point is located in the inside of either edge, determine a crossing point. */ /* 5) Take in turn each edge (vs,v_val0-2),...,(vs,v_val0 + val1-4) and (vs,v0), call the ith edge ei determine the closest point on each edge incident to v_val0+val1-3: */ /* Firstly determine numbers of elements in old and new configurations */ int edge_collapsible = 1, /* vertex stars */ *vstar0 = star0 + 2 * val0, *vstar1 = star1 + 2 * val1, i,j, i_min = -1; static Vertex vs[1], /* local copy of the simplified vertex */ a_triangle[3], wp[2], alpha[SIMPLIFY_VALENCE_MAX1]; /* val0 + val1 -4 <= valence max <=> val0, val1 <= valence_max +1 because val0 >=3 and val1 >=3 */ static float dp[2]; float tol_pivot, potential_err_pivot, largest_err_pivot = 0., d_min, lambda1, lambda2; _dxfErrorTestDataOnSurfaceInit(simp_data); /* set the maximum tolerance for the star */ tol_pivot = simp_data->tol_volume[vstar0[0]]; for (i=1;itol_volume[vstar0[i]] < tol_pivot) tol_pivot = simp_data->tol_volume[vstar0[i]]; for (i=0;itol_volume[vstar1[i]] < tol_pivot) tol_pivot = simp_data->tol_volume[vstar1[i]]; /* verify that the maximum tolerance is not already exceeded at the vertices */ i = 0; while (edge_collapsible && i < val0) /* the equality (star_err = tol) is allowed */ if (simp_data->err_volume[vstar0[i++]] > tol_pivot) edge_collapsible = 0; i = 0; while (edge_collapsible && i < val1) if (simp_data->err_volume[vstar1[i++]] > tol_pivot) edge_collapsible = 0; if (edge_collapsible) { Vertex *qi, *t0, *t1; /* local copy of the simplified vertex, translated by the due amount */ memcpy(vs[0], simp_data->simplified_vertex, sizeof (Vertex)); AddVec(vs[0], star_vert[val0+val1-2]); /* 1) */ /* then fabricate the triangle t0 = v1 v0 vval0-1 */ memcpy(a_triangle[0], simp_data->vert[vstar0[0]], sizeof (Vertex)); t0 = ExtractTriangleFromVertexStar(vstar1,0); /* find the point closest to vs on t0 and record the projection as well as the distance */ _dxfClosestPointOnTriangle(t0, vs[0], wp[0], alpha[0], &dp[0]); dp[0] += alpha[0][0] * simp_data->err_volume[vstar0[0]] + alpha[0][1] * simp_data->err_volume[vstar1[0]] + alpha[0][2] * simp_data->err_volume[vstar1[1]]; /* find the point closest to vs on t1 and record the same info triangle t1 = (v0, v1, v2) */ memcpy(a_triangle[0], simp_data->vert[vstar1[0]], sizeof (Vertex)); t1 = ExtractTriangleFromVertexStar(vstar0,0); _dxfClosestPointOnTriangle(t1, vs[0], wp[1], alpha[1], &dp[1]); dp[1] += alpha[1][0] * simp_data->err_volume[vstar1[0]] + alpha[1][1] * simp_data->err_volume[vstar0[0]] + alpha[1][2] * simp_data->err_volume[vstar0[1]]; /* the triangle with the smallest potential error wins: */ if (dp[0] < dp[1]) { /* use the barycentric coordinates to interpolate the potential data values */ _dxfErrorTestDataOnSurfaceInterpolate3(simp_data, vstar0[0], vstar1[0], vstar1[1], alpha[0][0], alpha[0][1], alpha[0][2]); potential_err_pivot = dp[0]; /* no need to test the data since we've just use this map point to compute new data values, so they have to match ! */ /* if (simp_data->vx_data) edge_collapsible = *//* SimpData simp_data, number_new_vertices_in_map, number_old_vertices, new_v2, new_v3, alpha_simp_vertex, alpha_n_v2, alpha_n_v3, old_v1, old_v2, old_v3, alpha_o_v1, alpha_o_v2, alpha_o_v3 *//* _dxfErrorTestDataOnSurfaceUpdate(simp_data, 1, 3, 0, 0, 1., 0., 0., vstar0[0], vstar1[0], vstar1[1], alpha[0][0], alpha[0][1], alpha[0][2]); */ } else { /* use the barycentric coordinates to interpolate the potential data values */ _dxfErrorTestDataOnSurfaceInterpolate3(simp_data, vstar1[0], vstar0[0], vstar0[1], alpha[1][0], alpha[1][1], alpha[1][2]); potential_err_pivot = dp[1]; /* no need to test the data since we've just use this map point to compute new data values, so they have to match ! */ } UPDATE_ERR_PIVOT; /* 2) determine the closest point to the vertex v0 on the triangles q1,...,qval0-2. */ if (edge_collapsible) { /* the triangles q1,...,val0-2 are constructed explicitely with the triples of vertices: (vs,v0,v1), vs(v_val03, v_val02) */ memcpy(a_triangle[0], vs[0], sizeof(Vertex)); i = 1; d_min = 0; while (i< val0-1) { qi = ExtractTriangleFromVertexStar(vstar0,i); _dxfClosestPointOnTriangle(qi, simp_data->vert[vstar1[0]], wp[0], alpha[i], &dp[0]); if (alpha[i][0]>0.) /* If v0 does not project on the link, then the projection is useful and can produce a valid constraint */ { dp[0] = ( dp[0] + simp_data->err_volume[vstar1[0]] - alpha[i][1] * simp_data->err_volume[vstar0[i]] - alpha[i][2] * simp_data->err_volume[vstar0[i+1]]) / alpha[i][0] ; if (i_min == -1 || dp[0] < d_min) {d_min = dp[0]; i_min = i;} } else { /* otherwise v0 projects onto the link edge of "triangle" the epsilon value at the simplified vertex cannot compensate but it is possible that the error values at the link are high enough to treat that constraint */ if (alpha[i][1] * simp_data->err_volume[vstar0[i]] + alpha[i][2] * simp_data->err_volume[vstar0[i+1]] >= dp[0]) { /* look no further for the closest triangle */ i_min = i; i = val0-1; /* next i will be = val0 and the loop will immediately stop */ d_min = 0.; } } i++; } /* the triangle corresponding to the smallest error wins */ if (i_min >=0) { /* if vertex0 does not project on the rim, i.e., if (alpha[i_min][0] > 0.) a potential new error value is generated at the pivot as follows: */ potential_err_pivot = d_min; UPDATE_ERR_PIVOT; if (simp_data->vx_data && edge_collapsible) edge_collapsible = _dxfErrorTestDataOnSurfaceUpdate(simp_data, 3, 1, vstar0[i_min], vstar0[i_min+1], alpha[i_min][0], alpha[i_min][1], alpha[i_min][2], vstar1[0], 0, 0, 1., 0., 0.); /* even though v0 == vstar1[0] is the same vertex, the data is allowed to be modified at that vertex, we still have to perform this test in case the simplified vertex would have been allowed to move: we would then verify that the old data is still compatible with the new directions */ } else edge_collapsible = 0; /* no projection was useful */ /* 4) */ if (edge_collapsible) { /* take each edge (vs,v1) .. (vs,v_val0-1) whereby v1, v_val0-1 refer to indices to the star the first vertex vstar0 */ i = 1; while (edge_collapsible && ivert[vstar0[i]], simp_data->vert[vstar1[0]], simp_data->vert[vstar0[j]], &lambda1, &lambda2, !j, 0, 1); /* If the closest point is located in the inside of either edge, determine a crossing point. */ if (lambda2 >= 0. && lambda2 <= 1.) /* inside 2 */ { d_min += lambda2 * simp_data->err_volume[vstar1[0]] + (1. -lambda2) * simp_data->err_volume[vstar0[j]]; if (simp_data->vx_data && lambda1 >= 0. && lambda1 <= 1.) edge_collapsible = _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 2, vstar0[i], 0, lambda1, 1. -lambda1, 0., vstar1[0], vstar0[j], 0, lambda2, 1. -lambda2, 0.); if (lambda1 > 0. && lambda1 <= 1.) /* inside 1 */ { potential_err_pivot = (d_min + ( lambda1 - 1.) * simp_data->err_volume[vstar0[i]]) / lambda1; UPDATE_ERR_PIVOT; } /* LAMBDA1 = 0 IS AGAIN A SPECIAL CASE if lambda1 = 0 we can't treat that constraint and the edge is not collapsible */ else if (lambda1 == 0.) { /* the error at vs cannot accomodate that constraint, but maybe it is already satisfied */ if (simp_data->err_volume[vstar0[i]] < d_min) edge_collapsible = 0; } } j++; if (i==j) j++; } i++; } /* 3) determine the closest point to the vertex v1 on the triangles qval0,...,q_val0+val1-3 */ if (edge_collapsible) { i_min = -1; /* the triangles qval0,...,q_val0+val1-3 are constructed explicitely using the vertices of vstar1 and the triples of vertices with the triples of vertices: (vs,v1,v2), (vs,v_val1-2, v_val1-1) */ i = 1; while (ivert[vstar0[0]], wp[0], alpha[i], &dp[0]); if (alpha[i][0]>0.) { dp[0] = ( dp[0] + simp_data->err_volume[vstar0[0]] - alpha[i][1] * simp_data->err_volume[vstar1[i]] - alpha[i][2] * simp_data->err_volume[vstar1[i+1]]) / alpha[i][0]; if (i_min == -1 || dp[0] < d_min) {d_min = dp[0]; i_min = i;} } else { /* otherwise v1 projects onto the link edge of "triangle" the epsilon value at the simplified vertex cannot compensate but it is possible that the error values at the link are high enough to treat that constraint */ if (alpha[i][1] * simp_data->err_volume[vstar1[i]] + alpha[i][2] * simp_data->err_volume[vstar1[i+1]] >= dp[0]) { /* look no further for the closest triangle */ i_min = i; i = val1-1; /* and the loop will immediately stop */ d_min = 0.; } } i++; } /* the closest triangle wins */ if (i_min >=0) { potential_err_pivot = d_min; UPDATE_ERR_PIVOT; if (simp_data->vx_data && edge_collapsible) edge_collapsible = _dxfErrorTestDataOnSurfaceUpdate(simp_data, 3, 1, vstar1[i_min], vstar1[i_min+1], alpha[i_min][0], alpha[i_min][1], alpha[i_min][2], vstar0[0], 0, 0, 1., 0., 0.); } else edge_collapsible = 0; /* no projection was useful */ /* 5) */ if (edge_collapsible) { /* Take in turn each edge (vs,v2),(vs,v1),..., (vs,v_val1-2) wherein v2... refer to star 1 original vertices of the star: 0,1,2,...,val1-1 we take 2,...,val1-2 */ i = 2; while (edge_collapsible && ivert[vstar1[i]], simp_data->vert[vstar0[0]], simp_data->vert[vstar1[j]], &lambda1, &lambda2, !j, 0, 1); if (lambda2 >= 0. && lambda2 <= 1.) /* inside 2 */ { d_min += lambda2 * simp_data->err_volume[vstar0[0]] + (1.-lambda2) * simp_data->err_volume[vstar1[j]]; if (simp_data->vx_data && lambda1 >= 0. && lambda1 <= 1.) edge_collapsible = _dxfErrorTestDataOnSurfaceUpdate(simp_data, 2, 2, vstar1[i], 0, lambda1, 1. -lambda1, 0., vstar0[0], vstar1[j], 0, lambda2, 1. -lambda2, 0.); if (lambda1 > 0. && lambda1 <= 1.) /* inside 1 */ { potential_err_pivot = (d_min - (1.-lambda1) * simp_data->err_volume[vstar1[i]]) / lambda1; UPDATE_ERR_PIVOT; } else if (lambda1 == 0.) { /* the error at vs cannot accomodate that constraint, but maybe it is already satisfied */ if (simp_data->err_volume[vstar1[i]] < d_min) edge_collapsible = 0; } } j++;if (i==j) j++; } i++; } } } } } } if (edge_collapsible) { /* update the tolerance and error values */ /* only the error at the pivot will change */ simp_data->err_volume[vstar1[0]] = largest_err_pivot; /* update the tolerance volume */ simp_data->tol_volume[vstar1[0]] = MIN(simp_data->tol_volume[vstar1[0]], simp_data->tol_volume[vstar0[0]]); _dxfErrorTestDataOnSurfaceEnd(simp_data, vstar1[0]); } return edge_collapsible; } /*+--------------------------------------------------------------------------+ | | | _dxfSimplifiedStarNormals | | | +--------------------------------------------------------------------------+*/ void _dxfSimplifiedStarNormals(Face *sface, Vertex *svert, Vertex *snor, float *sarea, float *scomp, Vertex simpvert, u_short first_face, u_short last_face) /*...........................................................................*/ /* compute the face normals in the simplified star as well as the compactness and area of the faces. */ /*...........................................................................*/ { register u_short i,i1; register u_char j,k; static VertexD triVert[3]; static double n[4]; /* introduce the simplified vertex */ CopyType2Type(float, double, simpvert, triVert[0], 3); for(i=first_face, i1=0; i < last_face; i++, i1++) { double sum_lengths_sq = 0.0; /* get the second and third vertices of the face.. */ u_short v1 = sface[i][1], v2 = sface[i][2]; /* ..and copy them at the 2nd and 3rd positions of faceVert */ CopyType2Type(float, double, svert[v1], triVert[1], 3); CopyType2Type(float, double, svert[v2], triVert[2], 3); _dxfTriangleNormalQR2D(triVert, n); CopyType2Type(double, float, n, snor[i1], 3); sarea[i1] = n[3]; /* compute the squared lenghts of the three edges: */ for(j=0;j<3;j++) { double len_edg_sq = 0, tmp; for(k=0;k<3;k++) { tmp = triVert[(j+1)%3][k] - triVert[j][k]; len_edg_sq += tmp * tmp; } sum_lengths_sq += len_edg_sq; } scomp[i1] = FOUR_SQRT_3 * sarea[i1] / sum_lengths_sq; } return; } /*+--------------------------------------------------------------------------+ | | | _dxfReplaceFaceNormals | | | +--------------------------------------------------------------------------+*/ void _dxfReplaceFaceNormals(int *star, Vertex *new_nor, float *new_area, float *new_comp, Vertex *nor, float *area, float *comp, u_short val) { /*...........................................................................*/ /* NOTE the count of the star goes from 0 to val-1 the corresponding new normals go from 2 to val-1 (triangle 0 and 1 are omitted since they will be removed) */ /*...........................................................................*/ register u_short i,i1; for(i=2, i1=0; itriangle_father[star1[0]] = star1[val1-1]; simp_data->triangle_father[star1[1]] = star1[2]; simp_data->vertex_father[v1] = v0; /* the remaining edges must be valid, and the way we choose to pick the triangle fathers dictates the way we choose edge father */ simp_data->edge_father[estar0[0]] = simp_data->edge_father[ estar1[val1-1] ] = estar0[1]; simp_data->edge_father[estar1[1]] = estar0[val0-1]; /* 2) change the valences of the vertices */ simp_data->valence[v0] = val0 + val1 - 4; /* v0 is not the only vertex whose valence changes: the valence of vstar0[1] and vstar1[1] decreases by one /\vstar0[1] | / \ | /____\ ---> | \ / | \ / | \/vstar1[1] | */ simp_data->valence[vstar0[1]] --; simp_data->valence[vstar1[1]] --; /* 3 remove the old edge lengths from the heap and put back the new ones */ /* 3.1 remove all edges in the star of v0 but v0 v1, because it was already removed from the heap (current edge under simplification) */ _dxfRemoveEdgesFromHeap(estar0+1, val0-1, simp_data->edge2index, simp_data); /* 3.2 remove all edges but v0-v1 in the star of v1 */ _dxfRemoveEdgesFromHeap(estar1+1, val1-1, simp_data->edge2index, simp_data); /* 4. mark the three edges that are to be discarded from the structure as collapsed with a -4 index */ /* update the number of edges left to be tested: the three collapsed edges wont be tested any more one among them was removed already : estar0[1] for the last two, we want to verify that they were not already counted as outside the heap */ if (simp_data->edge2index[estar1[val1-1]] != simp_data->heap_outside_index) simp_data->num_edg_remaining_4_testing --; if (simp_data->edge2index[estar1[1]] != simp_data->heap_outside_index) simp_data->num_edg_remaining_4_testing --; simp_data->edge2index[estar1[val1-1]] = simp_data->edge2index[estar1[1]] = -4; simp_data->edge2index[estar0[0]] = -7; /* 5.3 reinstate in the heap all the edges but 0 in estar0 that are their own parent, that is, all the edges that are a root of the disjoint set forest) */ n_edges_reinstated = _dxfReinstateEdgesInHeap(estar0+1, val0-1, simp_data, n_edges_added); /* 2.3 reinstate in the heap all edges but 0 and 1 and val1-1 in estar1 */ n_edges_reinstated += _dxfReinstateEdgesInHeap(estar1+2, val1-3, simp_data, n_edges_added+1); /* update the number of edges left to be tested */ simp_data->num_edg_remaining_4_testing += n_edges_reinstated; /* 3. copy the simplified vertex in the right place */ memcpy(simp_data->vert[v0], simp_data->simplified_vertex, sizeof(Vertex)); /* 4. replace the face normals, areas and compactness for the simplified stars with the new ones */ _dxfReplaceFaceNormals(star0, nor0, area0, comp0, simp_data->normal, simp_data->area, simp_data->compactness, val0); _dxfReplaceFaceNormals(star1, nor1, area1, comp1, simp_data->normal, simp_data->area, simp_data->compactness, val1); return n_edges_added[0] + n_edges_added[1]; } /*+--------------------------------------------------------------------------+ | | | _dxfSimplifyManifoldSurface | | | +--------------------------------------------------------------------------+*/ int _dxfSimplifyManifoldSurface(SimpData *simp_data) { int e, num_edg_tested = 0, /* number of edges extracted from the heap */ num_edg_added_to_heap = 0, num_edg_adjacent_2_boundary = 0; /* allocate buffers for storing various information on the vertex and edge stars */ /* since the valence of the simplified vertex can be as low as v0 + v1 - 4 (for an interior vertex) and is less than SIMPLIFY_VALENCE_MAX, assuming that v0 and v1 >= 1 we can be sure that the valence of v0 and v1 that are accepted for processing is <= SIMPLIFY_VALENCE_MAX+3. conservatively, we use SIMPLIFY_VALENCE_MAX+4 as a bound */ int vertex_star_buffer_size = 6 * SIMPLIFY_VALENCE_MAX4 * sizeof (int), edge_star_buffer_flsize = 9 * SIMPLIFY_VALENCE_MAX4 * sizeof (float), edge_star_buffer_vxsize = 4 * SIMPLIFY_VALENCE_MAX4 * sizeof (Vertex), edge_star_buffer_plsize = 4 * SIMPLIFY_VALENCE_MAX4 * sizeof (Plane), edge_star_buffer_fasize = 2 * SIMPLIFY_VALENCE_MAX4 * sizeof (Face), *vertex_star_buffer = NULL; Pointer edge_star_buffer_fl = NULL, edge_star_buffer_vx = NULL, edge_star_buffer_pl = NULL, edge_star_buffer_fa = NULL; vertex_star_buffer = (int *) DXAllocateZero(vertex_star_buffer_size); if (!vertex_star_buffer) goto error; edge_star_buffer_fl = (Pointer) DXAllocateZero(edge_star_buffer_flsize); if (!edge_star_buffer_fl) goto error; edge_star_buffer_vx = (Pointer) DXAllocateZero(edge_star_buffer_vxsize); if (!edge_star_buffer_vx) goto error; edge_star_buffer_pl = (Pointer) DXAllocateZero(edge_star_buffer_plsize); if (!edge_star_buffer_pl) goto error; edge_star_buffer_fa = (Pointer) DXAllocateZero(edge_star_buffer_fasize); if (!edge_star_buffer_fa) goto error; /* reference such storage arrays in the simp_data structure */ simp_data->vertex_star_buffer = vertex_star_buffer; simp_data->edge_star_buffer_fl = edge_star_buffer_fl; simp_data->edge_star_buffer_vx = edge_star_buffer_vx; /* 1) mark the edges touching the surface boundary as not simplifiable: */ if (!simp_data -> simplify_boundary) num_edg_adjacent_2_boundary = _dxfMarkEdgesAdjacent2Boundary(simp_data); /* 2) fill a Heap with edges, keyed with the weight of the edge */ while (num_edg_added_to_heap = _dxfBuildEdgeHeap(simp_data)){ simp_data->num_edg_weights += num_edg_added_to_heap; do { e = simp_data->get_lowest_weight_edge(simp_data); if (e != -1) { /* variable set to 1 if the edge is effectively collapsed */ int collapsible = 0; EdgeS *edg = simp_data->edge_array+e; int *endpoints = SurfEdgeGetVertices(edg), v0 = FATHER(endpoints[0], simp_data->vertex_father), v1 = FATHER(endpoints[1], simp_data->vertex_father); u_short val0 = simp_data->valence[v0], val1 = simp_data->valence[v1]; num_edg_tested ++; /* mark the edge as removed from the heap: give an index such that HeapDelete would fail + 1*/ simp_data->edge2index[e] = simp_data->heap_outside_index; /* decrement the number of edges remaining for testing */ simp_data->num_edg_remaining_4_testing --; /* treat a boundary vertex of the second kind */ if (simp_data->simplify_boundary && simp_data->boundary_vert[v0] && simp_data->boundary_vert[v1]) collapsible = _dxfCollapsibilityTestsBoundaryEdge2ndKind(simp_data, e, v0, v1, val0, val1, simp_data->preserve_boundary_length); /* treat a boundary vertex of the fisrt kind */ else if (simp_data->simplify_boundary && (simp_data->boundary_vert[v0] || simp_data->boundary_vert[v1])) collapsible = _dxfCollapsibilityTestsBoundaryEdge1stKind(simp_data, e, v0, v1, val0, val1); /* treat an interior edge */ else { /* perform the edge valence test for an interior edge */ if ( (val0 + val1 - 4 >= 3) && (val0 + val1 - 4 < simp_data->valence_max)) { /* BUILD THE 2 VERTEX STARS */ int *star0 = (int *) vertex_star_buffer, *estar0 = star0 + val0, *vstar0 = estar0 + val0, *star1 = vstar0 + val0, *estar1 = star1 + val1, *vstar1 = estar1 + val1; bzero(star0, 3 * (val0 + val1) * sizeof (int)); /* 3) construct the edge star */ _dxfBuildParentEdgeStars(edg,v0,v1,val0,val1,star0,star1,simp_data); /* 3.1) test whether the simplification is topologically feasible */ if (_dxfCollapseTopologicallyFeasible(simp_data, vstar0, vstar1, val0, val1)) { /* 3.2) compute a simplified vertex, and test whether the simplification is geometrically possible, if not, go to END */ /* the number of vertices is given by val0 + val1 - 2 and the number of triangles is also val0 + val1 - 2 */ u_short sval = val0 + val1 -2, sval1 = sval + 1, val02 = val0-2, val12 = val1-2; /* BUILD THE EDGE STAR */ Vertex *svert = (Vertex *) edge_star_buffer_vx, *snor0 = svert + sval1, *snor1 = snor0 + val02; Face *sface = (Face *) edge_star_buffer_fa; Plane *splane = (Plane *) edge_star_buffer_pl; float *scomp0 = (float *) edge_star_buffer_fl, *sarea0 = scomp0 + val02, *scomp1 = sarea0 + val02, *sarea1 = scomp1 + val12, *scomp = sarea1 + val12, *sarea = scomp + sval, *wspace = sarea + sval, min_comp_after; bzero((char *)svert, edge_star_buffer_vxsize); bzero((char *)sface, edge_star_buffer_fasize); bzero((char *)splane, edge_star_buffer_plsize); bzero((char *)scomp0, edge_star_buffer_flsize); /* the new origin is placed at the centroid of the edge link */ _dxfBuildEdgeStar(svert, sface, splane, sarea, scomp, v0, v1, val0, val1, star0, star1, simp_data); _dxfSimplifiedVertexLocation(simp_data, val0, val1, sval, svert, sface, splane, sarea, wspace, simp_data->preserve_volume); /* to determine whether the collapse is geometrically possible, the new vertex normals are needed */ _dxfSimplifiedStarNormals(sface, svert, snor0, sarea0, scomp0, simp_data->simplified_vertex, 1, val0-1); _dxfSimplifiedStarNormals(sface, svert, snor1, sarea1, scomp1, simp_data->simplified_vertex, val0, sval); /* minimum triangle compactness after simplification :*/ min_comp_after = MIN( _dxfMinFloatArray( scomp0, val02 ), _dxfMinFloatArray( scomp1, val12 ) ); if (_dxfCollapseGeometricallyFeasible(svert, sface, splane, scomp, min_comp_after, snor0, snor1, val0, sval, simp_data->min_scalprod_nor, simp_data->compactness_ratio)) { /* 3.3) if collapsible, test whether the new errors are within the tolerance limits */ collapsible = _dxfErrorWithinToleranceV(simp_data, svert, sface, splane, star0, star1, val0, val1); if (collapsible) { /* 3.4) simplify (collapse the edge ) and replace old edges in the heap with new ones */ simp_data->num_edg_collapsed += 3; /* translate back to the (0,0,0) origin */ AddVec( simp_data->simplified_vertex, svert[sval]); simp_data->num_edg_weights += _dxfCollapseEdge(simp_data, v0, v1, star0, star1, val0, val1, snor0, snor1, sarea0, sarea1, scomp0, scomp1); } else { /* rejected for tolerance off limits the edge can be reinserted afterwards if its length changes*/ simp_data->edg_rejected_4_tolerance ++; } } else {/* geometrically infeasible */ simp_data->edg_rejected_4_geometry ++; } } else { /* topologically infeasible */ simp_data->edge2index[e] = -3; /* this condition will not be modified even if the surface is simplified around the edge */ simp_data->edg_rejected_4_topology ++; } } else { /* If the valence is too high or too low, reject also for topology, However, this condition might change */ simp_data->edg_rejected_4_topology ++; } } /* end treat interior vertex */ } /* end if (e != -1) */ } while(e != -1); } /* end while (_dxfBuildEdgeHeap) */ DXFree((Pointer)vertex_star_buffer); DXFree(edge_star_buffer_fl); DXFree(edge_star_buffer_vx); DXFree(edge_star_buffer_pl); DXFree(edge_star_buffer_fa); return 1; error: if (vertex_star_buffer) DXFree((Pointer)vertex_star_buffer); if (edge_star_buffer_fl) DXFree(edge_star_buffer_fl); if (edge_star_buffer_vx) DXFree(edge_star_buffer_vx); if (edge_star_buffer_pl) DXFree(edge_star_buffer_pl); if (edge_star_buffer_fa) DXFree(edge_star_buffer_fa); return 0; } /*****************************************************************************/