--- ray/src/hd/sm_del.c 1999/03/05 16:32:22 3.11 +++ ray/src/hd/sm_del.c 2003/04/23 00:52:34 3.14 @@ -1,9 +1,6 @@ -/* Copyright (c) 1998 Silicon Graphics, Inc. */ - #ifndef lint -static char SCCSid[] = "$SunId$ SGI"; +static const char RCSid[] = "$Id: sm_del.c,v 3.14 2003/04/23 00:52:34 greg Exp $"; #endif - /* * sm_del.c */ @@ -11,34 +8,33 @@ static char SCCSid[] = "$SunId$ SGI"; #include "sm_flag.h" #include "sm_list.h" #include "sm_geom.h" -#include "sm_qtree.h" -#include "sm_stree.h" #include "sm.h" -static int Max_edges=200; +#define MAX_EDGE_INIT 100 +static int Max_edges= MAX_EDGE_INIT; static EDGE *Edges=NULL; static int Ecnt=0; -smFree_tri(sm,id) +smFree_tri(sm,id,t) SM *sm; int id; +TRI *t; { - TRI *tri,*t; - - tri = SM_NTH_TRI(sm,id); /* Add to the free_list */ - smClear_tri_flags(sm,id); - T_NEXT_FREE(tri) = SM_FREE_TRIS(sm); + T_NEXT_FREE(t) = SM_FREE_TRIS(sm); SM_FREE_TRIS(sm) = id; - T_VALID_FLAG(tri) = -1; +#ifdef DEBUG + T_VALID_FLAG(t) = INVALID; +#endif } /* Assumes mesh pointers have been cleaned up appropriately: just deletes from Point location and triangle data structure */ -smDelete_tri(sm,t_id) +smDelete_tri(sm,t_id,t) SM *sm; int t_id; +TRI *t; { /* NOTE: Assumes that a new triangle adjacent to each vertex @@ -47,16 +43,11 @@ int t_id; to id because the vertices can no longer point to tri id as being the first triangle pointer */ - SM_SAMPLE_TRIS(sm)--; - if(SM_IS_NTH_T_NEW(sm,t_id)) - smNew_tri_cnt--; + SM_CLR_NTH_T_ACTIVE(sm,t_id); + smFree_tri(sm,t_id,t); - smClear_tri_flags(sm,t_id); - - smFree_tri(sm,t_id); } - int eNew_edge() { @@ -66,35 +57,37 @@ eNew_edge() if(Ecnt >= Max_edges) { +#ifdef DEBUG if(Max_edges > 10000) - { - eputs("Too many edges in vertex loop\n"); - return(-1); - } + error(CONSISTENCY,"Too many edges in vertex loop\n"); +#endif Max_edges += 100; - if(!(Edges = (EDGE *)realloc(Edges,(Max_edges+1)*sizeof(EDGE)))) + if(!(Edges = (EDGE *)realloc((void *)Edges,(Max_edges+1)*sizeof(EDGE)))) goto memerr; } +#ifdef DEBUG + SET_E_NTH_TRI(Ecnt+1,0,INVALID); + SET_E_NTH_TRI(Ecnt+1,1,INVALID); +#endif return(++Ecnt); memerr: - error(SYSTEM,"eNew_edge(): Unable to allocate memory"); + error(SYSTEM,"eNew_edge): Unable to allocate memory"); } /* Return list of edges defining polygon formed by boundary of triangles adjacent to id. Return set of triangles adjacent to id to delete in delptr */ LIST -*smVertexStar(sm,id,del_ptr) +*smVertexStar(sm,id) SM *sm; -int id; -LIST **del_ptr; +S_ID id; { TRI *tri,*t_next; LIST *elist,*end; - int e,t_id,v_next,t_next_id,b_id,v_id; + int e,t_id,t_next_id,b_id,v_id,t_last_id; + S_ID v_next; - eClear_edges(); elist = end = NULL; /* Get the first triangle adjacent to vertex id */ @@ -102,7 +95,6 @@ LIST **del_ptr; tri = SM_NTH_TRI(sm,t_id); e = eNew_edge(); - /* Get the next vertex on the polygon boundary */ v_id = T_WHICH_V(tri,id); b_id = (v_id + 1)%3; @@ -110,21 +102,25 @@ LIST **del_ptr; SET_E_NTH_VERT(e,0,T_NTH_V(tri,b_id)); SET_E_NTH_TRI(e,0,INVALID); SET_E_NTH_TRI(e,1,T_NTH_NBR(tri,v_id)); - v_next = T_NTH_V(tri,(b_id+1)%3); + v_next = T_NTH_V(tri,(b_id+1)%3); SET_E_NTH_VERT(e,1,v_next); + elist = add_data_to_circular_list(elist,&end,e); t_next_id = t_id; t_next = tri; - - *del_ptr = push_data(*del_ptr,t_id); + t_last_id = -1; + /* Create a set to hold all of the triangles for deletion later */ while((t_next_id = T_NTH_NBR(t_next,b_id)) != t_id) { e = eNew_edge(); - if(e== INVALID) - return(NULL); + if(t_last_id != -1) + { + smDelete_tri(sm,t_last_id,t_next); + } t_next = SM_NTH_TRI(sm,t_next_id); + t_last_id = t_next_id; SET_E_NTH_VERT(e,0,v_next); SET_E_NTH_TRI(e,0,INVALID); v_id = T_WHICH_V(t_next,id); @@ -133,25 +129,23 @@ LIST **del_ptr; v_next = T_NTH_V(t_next,(b_id+1)%3); SET_E_NTH_VERT(e,1,v_next); elist = add_data_to_circular_list(elist,&end,e); - *del_ptr = push_data(*del_ptr,t_next_id); + } - return(elist); + smDelete_tri(sm,t_last_id,t_next); + smDelete_tri(sm,t_id,tri); + return(elist); } - int smTriangulate_add_tri(sm,id0,id1,id2,e0,e1,e2ptr) SM *sm; -int id0,id1,id2,e0,e1,*e2ptr; +S_ID id0,id1,id2; +int e0,e1,*e2ptr; { - int t_id; - int e2; + int t_id,e2; + TRI *t; -#ifdef DEBUG - if(id0 == INVALID || id1==INVALID || id2==INVALID) - error(CONSISTENCY,"bad id- smTriangulate_add_tri()\n"); -#endif - t_id = smAdd_tri(sm,id0,id1,id2); + t_id = smAdd_tri(sm,id0,id1,id2,&t); if(*e2ptr == 0) { e2 = eNew_edge(); @@ -164,48 +158,113 @@ int id0,id1,id2,e0,e1,*e2ptr; SET_E_NTH_TRI(e0,0,t_id); SET_E_NTH_TRI(e1,0,t_id); SET_E_NTH_TRI(e2,0,t_id); - +#ifdef DEBUG +#if DEBUG > 1 + if(E_NTH_TRI(e0,1) == E_NTH_TRI(e0,0) || + E_NTH_TRI(e1,1) == E_NTH_TRI(e1,0) || + E_NTH_TRI(e2,1) == E_NTH_TRI(e2,0)) + error(CONSISTENCY,"Non manifold\n"); +#endif +#endif *e2ptr = e2; return(t_id); + } -int -smTriangulateConvex(sm,plist,add_ptr) -SM *sm; -LIST *plist,**add_ptr; +int +smTriangulate_quad(sm,l) + SM *sm; + LIST *l; { - int t_id,e_id0,e_id1,e_id2; - int v_id0,v_id1,v_id2; - LIST *lptr; + int e1,e2,e3,e4,e,t_id; + S_ID id0,id1,id2,id3; + FVECT v0,v1,v2,v3,n; + double dp; + TRI *tri; + LIST *lptr; - lptr = plist; - e_id0 = (int)LIST_DATA(lptr); - v_id0 = E_NTH_VERT(e_id0,0); - lptr = LIST_NEXT(lptr); - while(LIST_NEXT(lptr) != plist) - { - e_id1 = (int)LIST_DATA(lptr); - v_id1 = E_NTH_VERT(e_id1,0); - v_id2 = E_NTH_VERT(e_id1,1); - lptr = LIST_NEXT(lptr); +#ifdef DEBUG + if(LIST_NEXT(LIST_NEXT(LIST_NEXT(LIST_NEXT(l)))) != l) + { + eputs("smTriangulate_quad(): not a quadrilateral\n"); + return(FALSE); + } + eputs("smTriangulate_quad():\n"); +#endif + lptr=l; + e1 = (int)LIST_DATA(l); + id0 = E_NTH_VERT(e1,0); + id1 = E_NTH_VERT(e1,1); + VSUB(v0,SM_NTH_WV(sm,id0),SM_VIEW_CENTER(sm)); + VSUB(v1,SM_NTH_WV(sm,id1),SM_VIEW_CENTER(sm)); + /* Get v2 */ + l = LIST_NEXT(l); + e2 = (int)LIST_DATA(l); + id2 = E_NTH_VERT(e2,1); + VSUB(v2,SM_NTH_WV(sm,id2),SM_VIEW_CENTER(sm)); + /* Get v3 */ + l = LIST_NEXT(l); + e3 = (int)LIST_DATA(l); + id3 = E_NTH_VERT(e3,1); + VSUB(v3,SM_NTH_WV(sm,id3),SM_VIEW_CENTER(sm)); + l = LIST_NEXT(l); + e4 = (int)LIST_DATA(l); + free_list(lptr); - if(LIST_NEXT(lptr) != plist) - e_id2 = 0; - else - e_id2 = (int)LIST_DATA(lptr); - t_id = smTriangulate_add_tri(sm,v_id0,v_id1,v_id2,e_id0,e_id1,&e_id2); - *add_ptr = push_data(*add_ptr,t_id); - e_id0 = -e_id2; + VCROSS(n,v0,v2); + normalize(n); + normalize(v1); + dp = DOT(n,v1); + e = 0; + if(dp > 0) + { + if(dp >= EV_EPS) + { + normalize(v3); + if(DOT(n,v3) < 0) + { + t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&e); + e *= -1; + t_id = smTriangulate_add_tri(sm,id2,id3,id0,e3,e4,&e); + return(TRUE); + } } - free_list(plist); - return(TRUE); -} -#ifdef TEST_DRIVER -FVECT Norm[500],B_V[500]; -int Ncnt,Bcnt,Del=0; + + } +#ifdef DEBUG +#if DEBUG > 1 + VCROSS(n,v1,v3); + normalize(n); + normalize(v0); + normalize(v2); + dp = DOT(n,v2); + if((dp < 0) || (dp < EV_EPS) || (DOT(n,v0) >= 0.0)) + error(CONSISTENCY,"smTriangulate_quad: cannot triangulate\n"); #endif +#endif + t_id = smTriangulate_add_tri(sm,id1,id2,id3,e2,e3,&e); + e *= -1; + t_id = smTriangulate_add_tri(sm,id3,id0,id1,e4,e1,&e); + return(TRUE); +} +eIn_tri(e,t) +int e; +TRI *t; +{ + if(T_NTH_V(t,0)==E_NTH_VERT(e,0)) + return(T_NTH_V(t,1)==E_NTH_VERT(e,1)||T_NTH_V(t,2)==E_NTH_VERT(e,1)); + else + if(T_NTH_V(t,1)==E_NTH_VERT(e,0)) + return(T_NTH_V(t,0)==E_NTH_VERT(e,1)||T_NTH_V(t,2)==E_NTH_VERT(e,1)); + else if(T_NTH_V(t,2)==E_NTH_VERT(e,0)) + return(T_NTH_V(t,0)==E_NTH_VERT(e,1)||T_NTH_V(t,1)==E_NTH_VERT(e,1)); + + return(FALSE); +} + + /* Triangulate the polygon defined by plist, and generating vertex p_id. Return list of added triangles in list add_ptr. Returns TRUE if successful, FALSE otherwise. This is NOT a general triangulation routine, @@ -213,36 +272,32 @@ int Ncnt,Bcnt,Del=0; */ int -smTriangulate(sm,id,plist,add_ptr) +smTriangulate(sm,id,plist) SM *sm; -int id; -LIST *plist,**add_ptr; +S_ID id; +LIST *plist; { LIST *l,*prev,*t; FVECT v0,v1,v2,n,p; - int is_tri,is_convex,cut,t_id,id0,id1,id2,e2,e1,enew; - double dp; - static int debug=0; + int is_tri,cut,t_id,e2,e1,enew; + S_ID id0,id1,id2; + double dp1,dp2; VSUB(p,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm)); - enew = 0; - is_convex = TRUE; - cut = is_tri= FALSE; - l = prev = plist; + normalize(p); + l = plist; + enew = 0; + cut = is_tri= FALSE; + prev = l; /* get v0,v1 */ e1 = (int)LIST_DATA(l); id0 = E_NTH_VERT(e1,0); id1 = E_NTH_VERT(e1,1); VSUB(v0,SM_NTH_WV(sm,id0),SM_VIEW_CENTER(sm)); + normalize(v0); VSUB(v1,SM_NTH_WV(sm,id1),SM_VIEW_CENTER(sm)); -#ifdef TEST_DRIVER - Del = TRUE; - VCOPY(B_V[0],v0); - VCOPY(B_V[1],v1); - Bcnt = 2; - Ncnt = 0; -#endif + normalize(v1); while(l) { l = LIST_NEXT(l); @@ -250,50 +305,48 @@ LIST *plist,**add_ptr; e2 = (int)LIST_DATA(l); id2 = E_NTH_VERT(e2,1); VSUB(v2,SM_NTH_WV(sm,id2),SM_VIEW_CENTER(sm)); -#ifdef TEST_DRIVER - VCOPY(B_V[Bcnt++],v2); -#endif + normalize(v2); if(LIST_NEXT(LIST_NEXT(l)) == prev) {/* Check if have a triangle */ is_tri = TRUE; break; } - /* determine if v0-v1-v2 is convex:defined clockwise on the sphere- so switch orientation */ - if(convex_angle(v2,v1,v0)) + VCROSS(n,v0,v2); + normalize(n); + dp1 = DOT(n,p); + if(dp1 <= 0.0) { /* test if safe to cut off v0-v1-v2 by testing if p lies outside of triangle v0-v1-v2: if so, because plist is the star polygon around p, the new edge v2-v0 cannot intersect any existing edges */ - VCROSS(n,v0,v2); - dp = DOT(n,p); - if(dp <= 0.0) - { + dp1 = DOT(n,v1); + /* Distance from point s to plane is d = |N.p0s|/||N|| */ + /* sp1 = s-p0 for point on plane p0, a.(b+c) = a.b + a.c */ + if(dp1 >= EV_EPS) + { /* remove edges e1,e2 and add triangle id0,id1,id2 */ - enew = 0; - t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew); - cut = TRUE; - *add_ptr = push_data(*add_ptr,t_id); - /* Insert edge enew into the list, reuse prev list element */ - LIST_NEXT(prev) = LIST_NEXT(l); - LIST_DATA(prev) = e1 = -enew; - /* If removing head of list- reset plist pointer */ - if(l== plist) - plist = prev; - /* free list element for e2 */ - LIST_NEXT(l)=NULL; - free_list(l); - l = prev; - VCOPY(v1,v2); - id1 = id2; - continue; + enew = 0; + t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew); + cut = TRUE; + /* Insert edge enew into the list, reuse prev list element */ + LIST_NEXT(prev) = LIST_NEXT(l); + LIST_DATA(prev) = e1 = -enew; + /* If removing head of list- reset plist pointer */ + if(l== plist) + plist = prev; + /* free list element for e2 */ + LIST_NEXT(l)=NULL; + free_list(l); + l = prev; + VCOPY(v1,v2); + id1 = id2; + continue; } - } - else - is_convex = FALSE; + } VCOPY(v0,v1); VCOPY(v1,v2); id0 = id1; @@ -308,40 +361,22 @@ LIST *plist,**add_ptr; is_tri = TRUE; break; } - - if(is_convex) - break; if(!cut) - { -#ifdef DEBUG - eputs("smTriangulate():Unable to triangulate\n"); -#endif - free_list(l); - while(*add_ptr) - { - t_id = pop_list(add_ptr); - smDelete_tri(sm,t_id); - } - return(FALSE); - } - + break; cut = FALSE; - is_convex = TRUE; - } - prev = l; } - if(is_tri) - { + prev = l; + } + if(is_tri) + { l = LIST_NEXT(l); enew = (int)LIST_DATA(l); t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew); - *add_ptr = push_data(*add_ptr,t_id); free_list(l); - } - else - if(!smTriangulateConvex(sm,l,add_ptr)) - return(FALSE); - + } + else + if(!smTriangulate_quad(sm,l)) + goto Terror; /* Set triangle adjacencies based on edge adjacencies */ FOR_ALL_EDGES(enew) { @@ -353,24 +388,242 @@ LIST *plist,**add_ptr; e2 = (T_WHICH_V(SM_NTH_TRI(sm,id1),E_NTH_VERT(enew,1))+2)%3; T_NTH_NBR(SM_NTH_TRI(sm,id1),e2) = id0; - } + } + +#ifdef DEBUG +#if DEBUG > 1 + { + TRI *t0,*t1,*n; + + FOR_ALL_EDGES(enew) + { + id0 = E_NTH_TRI(enew,0); + id1 = E_NTH_TRI(enew,1); + t0 = SM_NTH_TRI(sm,id0); + t1 = SM_NTH_TRI(sm,id1); + if(T_NTH_NBR(t0,0) == T_NTH_NBR(t0,1) || + T_NTH_NBR(t0,1) == T_NTH_NBR(t0,2) || + T_NTH_NBR(t0,0) == T_NTH_NBR(t0,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(t0,0)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(t0,1)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(t0,2)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); + + if(T_NTH_NBR(t1,0) == T_NTH_NBR(t1,1) || + T_NTH_NBR(t1,1) == T_NTH_NBR(t1,2) || + T_NTH_NBR(t1,0) == T_NTH_NBR(t1,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(t1,0)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(t1,1)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(t1,2)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); + } + + } +#endif +#endif + return(TRUE); + +Terror: + + error(CONSISTENCY,"smTriangulate():Unable to triangulate\n"); + } -eIn_tri(e,t) -int e; -TRI *t; + +void +smTris_swap_edge(sm,t_id,t1_id,e,e1,tn_id,tn1_id) + SM *sm; + int t_id,t1_id; + int e,e1; + int *tn_id,*tn1_id; { + S_ID verts[3]; + int enext,eprev,e1next,e1prev; + TRI *n,*ta,*tb,*t,*t1; + FVECT p1,p2,p3; + int ta_id,tb_id; + /* form new diagonal (e relative to t, and e1 relative to t1) + defined by quadrilateral formed by t,t1- swap for the opposite diagonal + */ + t = SM_NTH_TRI(sm,t_id); + t1 = SM_NTH_TRI(sm,t1_id); + enext = (e+1)%3; + eprev = (e+2)%3; + e1next = (e1+1)%3; + e1prev = (e1+2)%3; + verts[e] = T_NTH_V(t,e); + verts[enext] = T_NTH_V(t,enext); + verts[eprev] = T_NTH_V(t1,e1); + ta_id = smAdd_tri(sm,verts[0],verts[1],verts[2],&ta); +#if 0 + fprintf(stderr,"Added tri %d %d %d %d\n",ta_id,T_NTH_V(ta,0), + T_NTH_V(ta,1), T_NTH_V(ta,2)); +#endif + verts[e1] = T_NTH_V(t1,e1); + verts[e1next] = T_NTH_V(t1,e1next); + verts[e1prev] = T_NTH_V(t,e); + tb_id = smAdd_tri(sm,verts[0],verts[1],verts[2],&tb); +#if 0 + fprintf(stderr,"Added tri %d %d %d %d\n",tb_id,T_NTH_V(tb,0), + T_NTH_V(tb,1), T_NTH_V(tb,2)); +#endif + /* set the neighbors */ + T_NTH_NBR(ta,e) = T_NTH_NBR(t1,e1next); + T_NTH_NBR(tb,e1) = T_NTH_NBR(t,enext); + T_NTH_NBR(ta,enext)= tb_id; + T_NTH_NBR(tb,e1next)= ta_id; + T_NTH_NBR(ta,eprev)=T_NTH_NBR(t,eprev); + T_NTH_NBR(tb,e1prev)=T_NTH_NBR(t1,e1prev); - if(T_NTH_V(t,0)==E_NTH_VERT(e,0)) - return(T_NTH_V(t,1)==E_NTH_VERT(e,1)||T_NTH_V(t,2)==E_NTH_VERT(e,1)); - else - if(T_NTH_V(t,1)==E_NTH_VERT(e,0)) - return(T_NTH_V(t,0)==E_NTH_VERT(e,1)||T_NTH_V(t,2)==E_NTH_VERT(e,1)); - else if(T_NTH_V(t,2)==E_NTH_VERT(e,0)) - return(T_NTH_V(t,0)==E_NTH_VERT(e,1)||T_NTH_V(t,1)==E_NTH_VERT(e,1)); + /* Reset neighbor pointers of original neighbors */ + n = SM_NTH_TRI(sm,T_NTH_NBR(t,enext)); + T_NTH_NBR(n,T_NTH_NBR_PTR(t_id,n)) = tb_id; + n = SM_NTH_TRI(sm,T_NTH_NBR(t,eprev)); + T_NTH_NBR(n,T_NTH_NBR_PTR(t_id,n)) = ta_id; + n = SM_NTH_TRI(sm,T_NTH_NBR(t1,e1next)); + T_NTH_NBR(n,T_NTH_NBR_PTR(t1_id,n)) = ta_id; + n = SM_NTH_TRI(sm,T_NTH_NBR(t1,e1prev)); + T_NTH_NBR(n,T_NTH_NBR_PTR(t1_id,n)) = tb_id; - return(FALSE); + smDelete_tri(sm,t_id,t); + smDelete_tri(sm,t1_id,t1); + +#ifdef DEBUG +#if DEBUG > 1 + if(T_NTH_NBR(ta,0) == T_NTH_NBR(ta,1) || + T_NTH_NBR(ta,1) == T_NTH_NBR(ta,2) || + T_NTH_NBR(ta,0) == T_NTH_NBR(ta,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(ta,0)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(ta,1)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(ta,2)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); + if(T_NTH_NBR(ta,0) == T_NTH_NBR(ta,1) || + T_NTH_NBR(ta,1) == T_NTH_NBR(ta,2) || + T_NTH_NBR(ta,0) == T_NTH_NBR(ta,2)) + error(CONSISTENCY,"Invalid topology\n"); + + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(tb,0)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(tb,1)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); + n = SM_NTH_TRI(sm,T_NTH_NBR(tb,2)); + if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1) || + T_NTH_NBR(n,1) == T_NTH_NBR(n,2) || + T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) + error(CONSISTENCY,"Invalid topology\n"); + if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) || + !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2)))) + error(CONSISTENCY,"Invalid topology\n"); +#endif +#endif + *tn_id = ta_id; + *tn1_id = tb_id; + + return; } /* Test the new set of triangles for Delaunay condition. 'Edges' contains @@ -379,12 +632,12 @@ TRI *t; lies inside the circle defined by the CCW triangle- the edge is swapped for the opposite diagonal */ -smFixEdges(sm,add_list) +smFixEdges(sm) SM *sm; - LIST *add_list; { int e,t0_id,t1_id,e_new,e0,e1,e0_next,e1_next; - int i,v0_id,v1_id,v2_id,p_id,t0_nid,t1_nid; + S_ID v0_id,v1_id,v2_id,p_id; + int i,t0_nid,t1_nid; FVECT v0,v1,v2,p,np,v; TRI *t0,*t1; @@ -392,12 +645,10 @@ smFixEdges(sm,add_list) { t0_id = E_NTH_TRI(e,0); t1_id = E_NTH_TRI(e,1); - if((t0_id==INVALID) || (t1_id==INVALID)) - { #ifdef DEBUG - error(CONSISTENCY,"smFix_edges: Unassigned edge nbr\n"); + if((t0_id==INVALID) || (t1_id==INVALID)) + error(CONSISTENCY,"smFix_edges: Unassigned edge nbr\n"); #endif - } t0 = SM_NTH_TRI(sm,t0_id); t1 = SM_NTH_TRI(sm,t1_id); e0 = T_NTH_NBR_PTR(t1_id,t0); @@ -412,12 +663,11 @@ smFixEdges(sm,add_list) smDir(sm,v1,v1_id); smDir(sm,v2,v2_id); - VCOPY(p,SM_NTH_WV(sm,p_id)); - VSUB(p,p,SM_VIEW_CENTER(sm)); - if(point_in_cone(p,v0,v1,v2)) + VSUB(p,SM_NTH_WV(sm,p_id),SM_VIEW_CENTER(sm)); + normalize(p); + if(pt_in_cone(p,v2,v1,v0)) { - smTris_swap_edge(sm,t0_id,t1_id,e0,e1,&t0_nid,&t1_nid,&add_list); - + smTris_swap_edge(sm,t0_id,t1_id,e0,e1,&t0_nid,&t1_nid); /* Adjust the triangle pointers of the remaining edges to be processed */ @@ -438,6 +688,10 @@ smFixEdges(sm,add_list) else SET_E_NTH_TRI(i,1,t1_nid); } +#ifdef DEBUG + if(E_NTH_TRI(i,1) == E_NTH_TRI(i,0) ) + error(CONSISTENCY,"invalid edge\n"); +#endif } t0_id = t0_nid; t1_id = t1_nid; @@ -446,55 +700,64 @@ smFixEdges(sm,add_list) SET_E_NTH_VERT(e_new,1,v2_id); SET_E_NTH_TRI(e_new,0,t0_id); SET_E_NTH_TRI(e_new,1,t1_id); +#ifdef DEBUG + if(E_NTH_TRI(i,1) == E_NTH_TRI(i,0) ) + error(CONSISTENCY,"invalid edge\n"); +#endif } } - /* Add/Delete the appropriate triangles from the stree */ - smUpdate_locator(sm,add_list); } + +smDelete_samp(sm,s_id) +SM *sm; +S_ID s_id; +{ + QUADTREE qt; + S_ID *os; + + /* Mark as free */ + smUnalloc_samp(sm,s_id); + +#ifdef DEBUG + SM_NTH_VERT(sm,s_id) = INVALID; + /* fprintf(stderr,"deleting sample %d\n",s_id); */ +#endif + /* remove from its set */ + qt = SM_S_NTH_QT(sm,s_id); + os = qtqueryset(qt); + deletuelem(os, s_id); /* delete obj from unsorted os, no questions */ +} /* Remove vertex "id" from the mesh- and retriangulate the resulting hole: Returns TRUE if successful, FALSE otherwise. */ int smRemoveVertex(sm,id) SM *sm; - int id; + S_ID id; { - LIST *b_list,*add_list,*del_list; - int t_id,i; - static int cnt=0; - OBJECT *optr,*os; + LIST *b_list; /* generate list of edges that form the boundary of the - polygon formed by the triangles adjacent to vertex 'id' - */ - del_list = NULL; - b_list = smVertexStar(sm,id,&del_list); - if(!b_list) - { - if(del_list) - free_list(del_list); - return(FALSE); - } - add_list = NULL; + polygon formed by the triangles adjacent to vertex 'id'*/ + b_list = smVertexStar(sm,id); +#if 0 + {int i; + eputs("\n\n"); + for(i=1;i<=Ecnt;i++) + fprintf(stderr,"%d verts %d %d tris %d %d\n", + i,Edges[i].verts[0],Edges[i].verts[1], + Edges[i].tris[0],Edges[i].tris[1]); + } +#endif + /* Triangulate polygonal hole */ - if(!smTriangulate(sm,id,b_list,&add_list)) - { - free_list(del_list); - return(FALSE); - } - else - { - while(del_list) - { - t_id = pop_list(&del_list); - smDelete_tri(sm,t_id); - } - } - /* Fix up new triangles to be Delaunay-delnode contains set of - triangles to delete,add_list is the set of new triangles to add - */ - smFixEdges(sm,add_list); + smTriangulate(sm,id,b_list); + + /* Fix up new triangles to be Delaunay*/ + smFixEdges(sm); + smDelete_samp(sm,id); + eClear_edges(); return(TRUE); }