ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/sm_del.c
(Generate patch)

Comparing ray/src/hd/sm_del.c (file contents):
Revision 3.3 by gwlarson, Fri Sep 11 11:52:25 1998 UTC vs.
Revision 3.14 by greg, Wed Apr 23 00:52:34 2003 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1998 Silicon Graphics, Inc. */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ SGI";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /*
5   *  sm_del.c
6   */
7   #include "standard.h"
8 <
8 > #include "sm_flag.h"
9   #include "sm_list.h"
10   #include "sm_geom.h"
11   #include "sm.h"
12  
13 < static EDGE Edges[MAX_EDGES];
13 > #define MAX_EDGE_INIT 100
14 > static int Max_edges= MAX_EDGE_INIT;
15 > static EDGE *Edges=NULL;
16   static int Ecnt=0;
17  
18 < int
20 < remove_tri(qtptr,fptr,t_id)
21 <   QUADTREE *qtptr;
22 <   int *fptr;
23 <   int t_id;
24 < {
25 <    OBJECT tset[QT_MAXSET+1];  
26 <    int n;
27 <
28 <   if(QT_IS_EMPTY(*qtptr))
29 <     return(FALSE);
30 <   /* remove id from set */
31 <      else
32 <      {
33 <         if(!qtinset(*qtptr,t_id))
34 <           return(FALSE);
35 <         n = QT_SET_CNT(qtqueryset(*qtptr))-1;
36 <         *qtptr = qtdelelem(*qtptr,t_id);
37 <          if(n  == 0)
38 <             (*fptr) |= QT_COMPRESS;
39 <        if(!QT_FLAG_FILL_TRI(*fptr))
40 <          (*fptr)++;
41 <      }
42 <    return(TRUE);
43 < }
44 <
45 < int
46 < remove_tri_compress(qtptr,q0,q1,q2,t0,t1,t2,n,arg,t_id)
47 < QUADTREE *qtptr;
48 < FVECT q0,q1,q2;
49 < FVECT t0,t1,t2;
50 < int n;
51 < int *arg;
52 < int t_id;
53 < {
54 <    int f = 0;
55 <    /* NOTE compress */
56 <    return(remove_tri(qtptr,&f,t_id));
57 < }
58 <
59 <
60 <
61 < int
62 < stDelete_tri(st,t_id,t0,t1,t2)
63 <   STREE *st;
64 <   int t_id;
65 <   FVECT t0,t1,t2;
66 < {
67 <    int f;
68 <    FVECT dir;
69 <    
70 <  /* First add all of the leaf cells lying on the triangle perimeter:
71 <     mark all cells seen on the way
72 <   */
73 <    ST_CLEAR_FLAGS(st);
74 <    f = 0;
75 <    VSUB(dir,t1,t0);
76 <    stTrace_edge(st,t0,dir,1.0,remove_tri,&f,t_id);
77 <    VSUB(dir,t2,t1);
78 <    stTrace_edge(st,t1,dir,1.0,remove_tri,&f,t_id);
79 <    VSUB(dir,t0,t2);
80 <    stTrace_edge(st,t2,dir,1.0,remove_tri,&f,t_id);
81 <    /* Now visit interior */
82 <    if(QT_FLAG_FILL_TRI(f) || QT_FLAG_UPDATE(f))
83 <       stVisit_tri_interior(st,t0,t1,t2,remove_tri_compress,&f,t_id);
84 < }
85 <
86 <
87 < smLocator_remove_tri(sm,t_id)
18 > smFree_tri(sm,id,t)
19   SM *sm;
89 – int t_id;
90 – {
91 –  STREE *st;
92 –  char found;
93 –  TRI *t;
94 –  FVECT v0,v1,v2;
95 –
96 –  st = SM_LOCATOR(sm);
97 –
98 –  t = SM_NTH_TRI(sm,t_id);
99 –
100 –  VSUB(v0,SM_T_NTH_WV(sm,t,0),SM_VIEW_CENTER(sm));
101 –  VSUB(v1,SM_T_NTH_WV(sm,t,1),SM_VIEW_CENTER(sm));
102 –  VSUB(v2,SM_T_NTH_WV(sm,t,2),SM_VIEW_CENTER(sm));
103 –  found = stUpdate_tri(st,t_id,v0,v1,v2,remove_tri,remove_tri_compress);
104 –  return(found);
105 – }
106 –
107 – smFree_tri(sm,id)
108 – SM *sm;
20   int id;
21 + TRI *t;
22   {
111 –  TRI *tri;
112 –
113 –  tri = SM_NTH_TRI(sm,id);
23    /* Add to the free_list */
24 <  T_NEXT_FREE(tri) = SM_FREE_TRIS(sm);
24 >  T_NEXT_FREE(t) = SM_FREE_TRIS(sm);
25    SM_FREE_TRIS(sm) = id;
26 <  T_VALID_FLAG(tri) = -1;
26 > #ifdef DEBUG
27 >  T_VALID_FLAG(t) = INVALID;
28 > #endif
29   }
30  
31   /* Assumes mesh pointers have been cleaned up appropriately: just deletes from
32     Point location and triangle data structure
33   */
34 < smDelete_tri(sm,t_id)
34 > smDelete_tri(sm,t_id,t)
35   SM *sm;
36   int t_id;
37 + TRI *t;
38   {
39  
128 –
40    /* NOTE: Assumes that a new triangle adjacent to each vertex
41       has been added- before the deletion: replacing
42       the old tri- and therefore dont need to dereference any pointers
43       to id because the vertices can no longer
44       point to tri id as being the first triangle pointer
45    */
46 <  if(!SM_IS_NTH_T_BASE(sm,t_id))
47 <  {
137 <    SM_NUM_TRIS(sm)--;
138 <    if(SM_IS_NTH_T_NEW(sm,t_id))
139 <      smNew_tri_cnt--;
140 <  }
141 <  smClear_tri_flags(sm,t_id);
46 >  SM_CLR_NTH_T_ACTIVE(sm,t_id);
47 >  smFree_tri(sm,t_id,t);
48  
143 –  smFree_tri(sm,t_id);
144 –
49   }
50  
51 + int
52 + eNew_edge()
53 + {
54 +  if(!Edges)
55 +    if(!(Edges = (EDGE *)realloc(NULL,(Max_edges+1)*sizeof(EDGE))))
56 +      goto memerr;
57  
58 +  if(Ecnt >= Max_edges)
59 +    {
60 + #ifdef DEBUG
61 +      if(Max_edges > 10000)
62 +        error(CONSISTENCY,"Too many edges in vertex loop\n");
63 + #endif
64 +      Max_edges += 100;
65 +      if(!(Edges = (EDGE *)realloc((void *)Edges,(Max_edges+1)*sizeof(EDGE))))
66 +        goto memerr;
67 +    }
68 + #ifdef DEBUG
69 +  SET_E_NTH_TRI(Ecnt+1,0,INVALID);
70 +  SET_E_NTH_TRI(Ecnt+1,1,INVALID);
71 + #endif
72 +  return(++Ecnt);
73  
74 < LIST
75 < *smVertex_star_polygon(sm,id)
74 > memerr:
75 >  error(SYSTEM,"eNew_edge): Unable to allocate memory");
76 > }
77 >
78 > /* Return list of edges defining polygon formed by boundary of triangles
79 > adjacent to id. Return set of triangles adjacent to id to delete in delptr
80 > */
81 > LIST
82 > *smVertexStar(sm,id)
83   SM *sm;
84 < int id;
84 > S_ID id;
85   {
86      TRI *tri,*t_next;
87 <    LIST *elist,*end,*tlist;
88 <    int t_id,v_next,t_next_id;
89 <    int e;
87 >    LIST *elist,*end;
88 >    int e,t_id,t_next_id,b_id,v_id,t_last_id;
89 >    S_ID v_next;
90  
91      elist = end =  NULL;
92 +
93      /* Get the first triangle adjacent to vertex id */
94      t_id = SM_NTH_VERT(sm,id);
95      tri = SM_NTH_TRI(sm,t_id);
96  
97 <    
98 <    if((e = eNew_edge()) == SM_INVALID)
99 <    {
100 < #ifdef DEBUG
101 <        eputs("smVertex_star_polygon():Too many edges\n");
102 < #endif
103 <        return(NULL);
104 <    }
97 >    e = eNew_edge();
98 >    /* Get the  next vertex on the polygon boundary */
99 >    v_id = T_WHICH_V(tri,id);
100 >    b_id = (v_id + 1)%3;
101 >    /* Create an edge */
102 >    SET_E_NTH_VERT(e,0,T_NTH_V(tri,b_id));
103 >    SET_E_NTH_TRI(e,0,INVALID);
104 >    SET_E_NTH_TRI(e,1,T_NTH_NBR(tri,v_id));
105 >     v_next = T_NTH_V(tri,(b_id+1)%3);
106 >    SET_E_NTH_VERT(e,1,v_next);
107 >  
108      elist = add_data_to_circular_list(elist,&end,e);
173 –    v_next = (T_WHICH_V(tri,id)+1)%3;
174 –    SET_E_NTH_VERT(e,0,T_NTH_V(tri,v_next));
175 –    SET_E_NTH_TRI(e,0,SM_INVALID);
176 –    SET_E_NTH_TRI(e,1,T_NTH_NBR(tri,v_next));
177 –    v_next = (T_WHICH_V(tri,id)+2)%3;
178 –    SET_E_NTH_VERT(e,1,T_NTH_V(tri,v_next));
179 –
180 –    
109      t_next_id = t_id;
110      t_next = tri;
111 +    t_last_id = -1;
112 +
113 +    /* Create a set to hold all of the triangles for deletion later */
114  
115 <    tlist = push_data(NULL,t_id);
115 >    while((t_next_id = T_NTH_NBR(t_next,b_id)) != t_id)
116 >    {
117 >      e = eNew_edge();
118 >      if(t_last_id != -1)
119 >      {
120 >        smDelete_tri(sm,t_last_id,t_next);
121 >      }
122 >      t_next = SM_NTH_TRI(sm,t_next_id);
123 >      t_last_id = t_next_id;
124 >      SET_E_NTH_VERT(e,0,v_next);
125 >      SET_E_NTH_TRI(e,0,INVALID);
126 >      v_id = T_WHICH_V(t_next,id);
127 >      b_id = (v_id + 1)%3;
128 >      SET_E_NTH_TRI(e,1,T_NTH_NBR(t_next,v_id));
129 >      v_next = T_NTH_V(t_next,(b_id+1)%3);
130 >      SET_E_NTH_VERT(e,1,v_next);
131 >      elist = add_data_to_circular_list(elist,&end,e);
132  
186 –    while((t_next_id = smTri_next_ccw_nbr(sm,t_next,id)) != t_id)
187 –    {  
188 –        if((e = eNew_edge()) == SM_INVALID)
189 –        {
190 – #ifdef DEBUG
191 –            eputs("smVertex_star_polygon():Too many edges\n");
192 – #endif
193 –            return(NULL);
194 –        }
195 –        elist = add_data_to_circular_list(elist,&end,e);
196 –        t_next = SM_NTH_TRI(sm,t_next_id);
197 –        v_next = (T_WHICH_V(t_next,id)+1)%3;
198 –        SET_E_NTH_VERT(e,0,T_NTH_V(t_next,v_next));
199 –        SET_E_NTH_TRI(e,0,SM_INVALID);
200 –        SET_E_NTH_TRI(e,1,T_NTH_NBR(t_next,v_next));
201 –        v_next = (T_WHICH_V(t_next,id)+2)%3;
202 –        SET_E_NTH_VERT(e,1,T_NTH_V(t_next,v_next));
203 –        tlist = push_data(tlist,t_next_id);
133      }
134 <    while(tlist)
135 <    {
136 <        t_id = (int)pop_list(&tlist);
208 <        /* first remove from point location structure */
209 <        smLocator_remove_tri(sm,t_id);
210 <        smDelete_tri(sm,t_id);
211 <    }
212 <    return(elist);
134 >    smDelete_tri(sm,t_last_id,t_next);
135 >    smDelete_tri(sm,t_id,tri);
136 >   return(elist);
137   }
138  
139   int
140 < smEdge_intersect_polygon(sm,v0,v1,l)
140 > smTriangulate_add_tri(sm,id0,id1,id2,e0,e1,e2ptr)
141   SM *sm;
142 < FVECT v0,v1;
143 < LIST *l;
142 > S_ID id0,id1,id2;
143 > int e0,e1,*e2ptr;
144   {
145 <    FVECT e0,e1;
146 <    int e,id_e0,id_e1;
223 <    LIST *el,*eptr;
224 <    
225 <    /* Test the edges in l against v0v1 to see if v0v1 intersects
226 <       any other edges
227 <     */
228 <    
229 <    el = l;
145 >  int t_id,e2;
146 >  TRI *t;
147  
148 <    while(el)
149 <    {
150 <      e = (int)LIST_DATA(el);
151 <      id_e0 = E_NTH_VERT(e,0);
152 <      id_e1 = E_NTH_VERT(e,1);
153 <      /* NOTE: DO these need to be normalized? Just subtract center? */
154 <      smDir(sm,e0,id_e0);
155 <      smDir(sm,e1,id_e1);
156 <      if(sedge_intersect(v0,v1,e0,e1))
157 <        return(TRUE);
148 >  t_id = smAdd_tri(sm,id0,id1,id2,&t);
149 >  if(*e2ptr == 0)
150 >  {
151 >    e2 = eNew_edge();
152 >    SET_E_NTH_VERT(e2,0,id2);
153 >    SET_E_NTH_VERT(e2,1,id0);
154 >  }
155 >  else
156 >    e2 = *e2ptr;
157 >  /* set appropriate tri for each edge*/
158 >  SET_E_NTH_TRI(e0,0,t_id);
159 >  SET_E_NTH_TRI(e1,0,t_id);
160 >  SET_E_NTH_TRI(e2,0,t_id);
161 > #ifdef DEBUG
162 > #if DEBUG > 1
163 >  if(E_NTH_TRI(e0,1) == E_NTH_TRI(e0,0) ||
164 >     E_NTH_TRI(e1,1) == E_NTH_TRI(e1,0) ||
165 >     E_NTH_TRI(e2,1) == E_NTH_TRI(e2,0))
166 >    error(CONSISTENCY,"Non manifold\n");
167 > #endif
168 > #endif
169 >  *e2ptr = e2;
170 >  return(t_id);
171  
242 –      el = LIST_NEXT(el);
243 –      if(el == l)
244 –        break;
245 –    }
246 –    return(FALSE);
172   }
173  
174 < int
175 < smFind_next_convex_vertex(sm,id0,id1,v0,v1,l)
176 <   SM *sm;
177 <   int id0,id1;
253 <   FVECT v0,v1;
254 <   LIST *l;
174 > int
175 > smTriangulate_quad(sm,l)
176 >     SM *sm;
177 >     LIST *l;
178   {
179 <    int e,id;
180 <    LIST *el;
181 <    FVECT v;
179 >  int e1,e2,e3,e4,e,t_id;
180 >  S_ID id0,id1,id2,id3;
181 >  FVECT v0,v1,v2,v3,n;
182 >  double dp;
183 >  TRI *tri;
184 >  LIST *lptr;
185  
260 –    /* starting with the end of edge at head of l, search sequentially for
261 –      vertex v such that v0v1v is a convex angle, and the edge v1v does
262 –      not intersect any other edges
263 –   */
264 –    id = SM_INVALID;
265 –    el = l;
266 –    while(id != id0)
267 –    {
268 –        e = (int)LIST_DATA(el);
269 –        id = E_NTH_VERT(e,1);
270 –
271 –        smDir(sm,v,id);
272 –
273 –        if(convex_angle(v0,v1,v) && !smEdge_intersect_polygon(sm,v1,v,l))
274 –           return(id);
275 –              
276 –        el = LIST_NEXT(el);
277 –        if(el == l)
278 –           break;
279 –    }
280 –    return(SM_INVALID);
281 – }
282 –
283 – int
284 – split_edge_list(id0,id_new,l,lnew)
285 – int id0,id_new;
286 – LIST **l,**lnew;
287 – {
288 –    LIST *list,*lptr,*end;
289 –    int e,e1,e2,new_e;
290 –
291 –    e2 = SM_INVALID;
292 –    list = lptr = *l;
293 –
294 –    if((new_e = eNew_edge())==SM_INVALID)
295 –     {
186   #ifdef DEBUG
187 <            eputs("split_edge_list():Too many edges\n");
187 >  if(LIST_NEXT(LIST_NEXT(LIST_NEXT(LIST_NEXT(l)))) != l)
188 >  {
189 >    eputs("smTriangulate_quad(): not a quadrilateral\n");
190 >    return(FALSE);
191 >  }
192 >    eputs("smTriangulate_quad():\n");
193   #endif
194 <         return(FALSE);
195 <     }
196 <    SET_E_NTH_VERT(new_e,0,id0);
197 <    SET_E_NTH_VERT(new_e,1,id_new);
198 <    SET_E_NTH_TRI(new_e,0,SM_INVALID);
199 <    SET_E_NTH_TRI(new_e,1,SM_INVALID);
200 <    
201 <    while(e2 != id_new)
202 <    {
203 <        lptr = LIST_NEXT(lptr);
204 <        e = (int)LIST_DATA(lptr);
205 <        e2 = E_NTH_VERT(e,1);
206 <        if(lptr == list)
207 <        {
208 < #ifdef DEBUG        
209 <          eputs("split_edge_list():cant find vertex\n");
210 < #endif
211 <          *lnew = NULL;
212 <          return(FALSE);
318 <        }
194 >  lptr=l;
195 >  e1 = (int)LIST_DATA(l);
196 >  id0 = E_NTH_VERT(e1,0);
197 >  id1 = E_NTH_VERT(e1,1);
198 >  VSUB(v0,SM_NTH_WV(sm,id0),SM_VIEW_CENTER(sm));
199 >  VSUB(v1,SM_NTH_WV(sm,id1),SM_VIEW_CENTER(sm));  
200 >  /* Get v2 */
201 >  l = LIST_NEXT(l);
202 >  e2 = (int)LIST_DATA(l);
203 >  id2 = E_NTH_VERT(e2,1);
204 >  VSUB(v2,SM_NTH_WV(sm,id2),SM_VIEW_CENTER(sm));
205 >  /* Get v3 */
206 >  l = LIST_NEXT(l);
207 >  e3 = (int)LIST_DATA(l);
208 >  id3 = E_NTH_VERT(e3,1);
209 >  VSUB(v3,SM_NTH_WV(sm,id3),SM_VIEW_CENTER(sm));
210 >  l = LIST_NEXT(l);
211 >  e4 = (int)LIST_DATA(l);
212 >  free_list(lptr);
213  
214 +  VCROSS(n,v0,v2);
215 +  normalize(n);
216 +  normalize(v1);
217 +  dp = DOT(n,v1);
218 +  e = 0;
219 +  if(dp > 0)
220 +  {
221 +    if(dp  >= EV_EPS)
222 +    {
223 +      normalize(v3);
224 +      if(DOT(n,v3) < 0)
225 +      {
226 +        t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&e);
227 +        e *= -1;
228 +        t_id = smTriangulate_add_tri(sm,id2,id3,id0,e3,e4,&e);
229 +        return(TRUE);
230 +      }
231      }
321 –    end = lptr;
322 –    lptr =  LIST_NEXT(lptr);
323 –    list = add_data_to_circular_list(list,&end,-new_e);
324 –    *lnew = list;
232  
233 <    /* now follow other cycle */
234 <
235 <    list = lptr;
236 <    e2 = SM_INVALID;
237 <    while(e2 != id0)
238 <    {
239 <        lptr = LIST_NEXT(lptr);
240 <        e = (int)LIST_DATA(lptr);
241 <        e2 = E_NTH_VERT(e,1);
242 <        if(lptr == list)
336 <        {
337 < #ifdef DEBUG        
338 <          eputs("split_edge_list():cant find intial vertex\n");
233 >  }
234 > #ifdef DEBUG
235 > #if DEBUG > 1
236 >  VCROSS(n,v1,v3);
237 >  normalize(n);
238 >  normalize(v0);
239 >  normalize(v2);
240 >  dp = DOT(n,v2);
241 >  if((dp < 0) || (dp < EV_EPS) || (DOT(n,v0) >= 0.0))  
242 >    error(CONSISTENCY,"smTriangulate_quad: cannot triangulate\n");
243   #endif
244 <          *l = NULL;
245 <          return(FALSE);
246 <        }
247 <
248 <    }
345 <    end = lptr;
346 <    list = add_data_to_circular_list(list,&end,new_e);
347 <    *l = list;
348 <    return(TRUE);
244 > #endif
245 >  t_id = smTriangulate_add_tri(sm,id1,id2,id3,e2,e3,&e);
246 >  e *= -1;
247 >  t_id = smTriangulate_add_tri(sm,id3,id0,id1,e4,e1,&e);
248 >  return(TRUE);
249   }
250  
251 <
252 < int
253 < smTriangulate_convex(sm,plist)
354 < SM *sm;
355 < LIST *plist;
251 > eIn_tri(e,t)
252 > int e;
253 > TRI *t;
254   {
357 –    TRI *tri;
358 –    int t_id,e_id0,e_id1,e_id2;
359 –    int v_id0,v_id1,v_id2;
360 –    LIST *lptr;
361 –    int cnt;
255  
256 <    lptr = plist;
257 <    e_id0 = (int)LIST_DATA(lptr);
258 <    v_id0 = E_NTH_VERT(e_id0,0);
259 <    lptr = LIST_NEXT(lptr);
260 <    while(LIST_NEXT(lptr) != plist)
261 <    {
262 <        e_id1 = (int)LIST_DATA(lptr);
370 <        v_id1 = E_NTH_VERT(e_id1,0);
371 <        v_id2 = E_NTH_VERT(e_id1,1);
372 <        /* form a triangle for each triple of with v0 as base of star */
373 <        t_id = smAdd_tri(sm,v_id0,v_id1,v_id2,&tri);
374 <        smLocator_add_tri(sm,t_id,v_id0,v_id1,v_id2);
375 <        /* add which pointer?*/
256 >  if(T_NTH_V(t,0)==E_NTH_VERT(e,0))
257 >    return(T_NTH_V(t,1)==E_NTH_VERT(e,1)||T_NTH_V(t,2)==E_NTH_VERT(e,1));
258 >  else
259 >    if(T_NTH_V(t,1)==E_NTH_VERT(e,0))
260 >      return(T_NTH_V(t,0)==E_NTH_VERT(e,1)||T_NTH_V(t,2)==E_NTH_VERT(e,1));
261 >    else if(T_NTH_V(t,2)==E_NTH_VERT(e,0))
262 >      return(T_NTH_V(t,0)==E_NTH_VERT(e,1)||T_NTH_V(t,1)==E_NTH_VERT(e,1));
263  
264 <        lptr = LIST_NEXT(lptr);
264 >  return(FALSE);
265 > }
266  
379 –        if(LIST_NEXT(lptr) != plist)
380 –        {
381 –            e_id2 = eNew_edge();
382 –            SET_E_NTH_VERT(e_id2,0,v_id2);
383 –            SET_E_NTH_VERT(e_id2,1,v_id0);
384 –        }
385 –        else
386 –           e_id2 = (int)LIST_DATA(lptr);
387 –        
388 –        /* set appropriate tri for each edge*/
389 –        SET_E_NTH_TRI(e_id0,0,t_id);
390 –        SET_E_NTH_TRI(e_id1,0,t_id);
391 –        SET_E_NTH_TRI(e_id2,0,t_id);
267  
268 <        e_id0 = -e_id2;
269 <    }
270 <    
271 <    free_list(plist);
272 <    return(TRUE);
273 < }
268 > /* Triangulate the polygon defined by plist, and generating vertex p_id.
269 >   Return list of added triangles in list add_ptr. Returns TRUE if
270 >   successful, FALSE otherwise. This is NOT a general triangulation routine,
271 >   assumes polygon star relative to id
272 > */
273 >
274   int
275 < smTriangulate_elist(sm,plist)
275 > smTriangulate(sm,id,plist)
276   SM *sm;
277 + S_ID id;
278   LIST *plist;
279   {
280 <    LIST *l,*el1;
281 <    FVECT v0,v1,v2;
282 <    int id0,id1,id2,e,id_next;
283 <    char flipped;
284 <    int done;
280 >    LIST *l,*prev,*t;
281 >    FVECT v0,v1,v2,n,p;
282 >    int is_tri,cut,t_id,e2,e1,enew;
283 >    S_ID id0,id1,id2;
284 >    double dp1,dp2;
285  
286 +    VSUB(p,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
287 +    normalize(p);
288      l = plist;
289 <    
289 >
290 >    enew = 0;
291 >    cut = is_tri=  FALSE;
292 >    prev = l;
293 >    /* get v0,v1 */
294 >    e1 = (int)LIST_DATA(l);
295 >    id0 = E_NTH_VERT(e1,0);
296 >    id1 = E_NTH_VERT(e1,1);
297 >    VSUB(v0,SM_NTH_WV(sm,id0),SM_VIEW_CENTER(sm));
298 >    normalize(v0);
299 >    VSUB(v1,SM_NTH_WV(sm,id1),SM_VIEW_CENTER(sm));  
300 >    normalize(v1);
301      while(l)
302      {
414 –        /* get v0,v1,v2 */
415 –      e = (int)LIST_DATA(l);
416 –      id0 = E_NTH_VERT(e,0);
417 –      id1 = E_NTH_VERT(e,1);
303        l = LIST_NEXT(l);
304 <      e = (int)LIST_DATA(l);
305 <      id2 = E_NTH_VERT(e,1);
306 <
307 <      smDir(sm,v0,id0);
308 <      smDir(sm,v1,id1);
309 <      smDir(sm,v2,id2);
310 <      /* determine if convex (left turn), or concave(right turn) angle */
311 <      if(convex_angle(v0,v1,v2))
312 <      {
428 <        if(l == plist)
429 <          break;
430 <        else
431 <          continue;
304 >      /* Get v2 */
305 >      e2 = (int)LIST_DATA(l);
306 >      id2 = E_NTH_VERT(e2,1);
307 >      VSUB(v2,SM_NTH_WV(sm,id2),SM_VIEW_CENTER(sm));
308 >      normalize(v2);
309 >      if(LIST_NEXT(LIST_NEXT(l)) == prev)
310 >      {/* Check if have a triangle */
311 >           is_tri = TRUE;      
312 >           break;
313        }
314 <      /* if concave: add edge and recurse on two sub polygons */
315 <      id_next = smFind_next_convex_vertex(sm,id0,id1,v0,v1,LIST_NEXT(l));
316 <      if(id_next == SM_INVALID)
314 >      /* determine if v0-v1-v2 is convex:defined clockwise on the sphere-
315 >       so switch orientation
316 >       */
317 >      VCROSS(n,v0,v2);
318 >      normalize(n);
319 >      dp1 = DOT(n,p);      
320 >      if(dp1  <= 0.0)
321        {
322 < #ifdef DEBUG
323 <          eputs("smTriangulate_elist():Unable to find convex vertex\n");
324 < #endif
325 <          return(FALSE);
326 <      }
327 <      /* add edge */
328 <      el1 = NULL;
329 <      /* Split edge list l into two lists: one from id1-id_next-id1,
330 <         and the next from id2-id_next-id2
331 <      */
332 <      split_edge_list(id1,id_next,&l,&el1);
333 <      /* Recurse and triangulate the two edge lists */
334 <      done = smTriangulate_elist(sm,l);
335 <      if(done)
336 <        done = smTriangulate_elist(sm,el1);
337 <      return(done);
322 >          /* test if safe to cut off v0-v1-v2 by testing if p lies outside of
323 >         triangle v0-v1-v2: if so, because plist is the star polygon around p,
324 >          the new edge v2-v0 cannot intersect any existing edges
325 >        */
326 >        dp1 = DOT(n,v1);
327 >        /* Distance from point s to plane is d = |N.p0s|/||N|| */
328 >        /* sp1 = s-p0 for point on plane p0, a.(b+c) = a.b + a.c */
329 >        if(dp1 >= EV_EPS)
330 >        {
331 >            /* remove edges e1,e2 and add triangle id0,id1,id2 */
332 >            enew = 0;
333 >            t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew);
334 >            cut = TRUE;
335 >            /* Insert edge enew into the list, reuse prev list element */
336 >            LIST_NEXT(prev) = LIST_NEXT(l);
337 >            LIST_DATA(prev) = e1 = -enew;
338 >            /* If removing head of list- reset plist pointer */
339 >            if(l== plist)
340 >              plist = prev;
341 >            /* free list element for e2 */
342 >            LIST_NEXT(l)=NULL;
343 >            free_list(l);
344 >            l = prev;
345 >            VCOPY(v1,v2);
346 >            id1 = id2;
347 >            continue;
348 >        }
349      }
350 <    done = smTriangulate_convex(sm,plist);
351 <    return(done);
352 < }
350 >      VCOPY(v0,v1);
351 >      VCOPY(v1,v2);
352 >      id0 = id1;
353 >      id1 = id2;
354 >      e1 = e2;  
355 >      /* check if gone around circular list without adding any
356 >         triangles: prevent infinite loop */
357 >      if(l == plist)
358 >      {
359 >        if(LIST_NEXT(LIST_NEXT(l)) == prev)
360 >          {/* Check if have a triangle */
361 >            is_tri = TRUE;      
362 >            break;
363 >          }
364 >        if(!cut)
365 >          break;
366 >        cut = FALSE;
367 >    }
368 >      prev = l;
369 >  }
370 >  if(is_tri)
371 >  {
372 >      l = LIST_NEXT(l);
373 >      enew = (int)LIST_DATA(l);
374 >      t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew);
375 >      free_list(l);
376 >  }
377 >  else
378 >     if(!smTriangulate_quad(sm,l))
379 >        goto Terror;
380 >    /* Set triangle adjacencies based on edge adjacencies */
381 >    FOR_ALL_EDGES(enew)
382 >    {
383 >      id0 = E_NTH_TRI(enew,0);
384 >      id1 = E_NTH_TRI(enew,1);
385 >        
386 >      e1 = (T_WHICH_V(SM_NTH_TRI(sm,id0),E_NTH_VERT(enew,0))+2)%3;
387 >      T_NTH_NBR(SM_NTH_TRI(sm,id0),e1) = id1;
388 >      
389 >      e2 = (T_WHICH_V(SM_NTH_TRI(sm,id1),E_NTH_VERT(enew,1))+2)%3;
390 >      T_NTH_NBR(SM_NTH_TRI(sm,id1),e2) = id0;
391 >  }
392  
393 < int
394 < smTriangulate(sm,plist)
395 < SM *sm;
396 < LIST *plist;
462 < {
463 <    int e,id_t0,id_t1,e0,e1;
464 <    TRI *t0,*t1;
465 <    int test;
466 <    
467 <    test = smTriangulate_elist(sm,plist);
393 > #ifdef DEBUG
394 > #if DEBUG > 1
395 >    {
396 >      TRI *t0,*t1,*n;
397  
398 <    if(!test)
470 <       return(test);
471 <    FOR_ALL_EDGES(e)
398 >    FOR_ALL_EDGES(enew)
399      {
400 <        id_t0 = E_NTH_TRI(e,0);
401 <        id_t1 = E_NTH_TRI(e,1);
402 <        if((id_t0==SM_INVALID) || (id_t1==SM_INVALID))
403 <        {
404 < #ifdef DEBUG
405 <           eputs("smTriangulate(): Unassigned edge neighbor\n");
400 >      id0 = E_NTH_TRI(enew,0);
401 >      id1 = E_NTH_TRI(enew,1);
402 >      t0 = SM_NTH_TRI(sm,id0);
403 >      t1 = SM_NTH_TRI(sm,id1);
404 >      if(T_NTH_NBR(t0,0) == T_NTH_NBR(t0,1)  ||
405 >       T_NTH_NBR(t0,1) == T_NTH_NBR(t0,2)  ||
406 >       T_NTH_NBR(t0,0) == T_NTH_NBR(t0,2))
407 >      error(CONSISTENCY,"Invalid topology\n");
408 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,0))) ||
409 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,1))) ||
410 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,2))))
411 >      error(CONSISTENCY,"Invalid topology\n");
412 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,0));
413 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
414 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
415 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
416 >      error(CONSISTENCY,"Invalid topology\n");
417 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
418 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
419 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
420 >      error(CONSISTENCY,"Invalid topology\n");
421 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,1));
422 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
423 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
424 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
425 >      error(CONSISTENCY,"Invalid topology\n");
426 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
427 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
428 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
429 >      error(CONSISTENCY,"Invalid topology\n");
430 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,2));
431 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
432 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
433 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
434 >      error(CONSISTENCY,"Invalid topology\n");
435 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
436 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
437 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
438 >      error(CONSISTENCY,"Invalid topology\n");
439 >
440 >      if(T_NTH_NBR(t1,0) == T_NTH_NBR(t1,1)  ||
441 >       T_NTH_NBR(t1,1) == T_NTH_NBR(t1,2)  ||
442 >       T_NTH_NBR(t1,0) == T_NTH_NBR(t1,2))
443 >      error(CONSISTENCY,"Invalid topology\n");
444 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,0))) ||
445 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,1))) ||
446 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,2))))
447 >      error(CONSISTENCY,"Invalid topology\n");
448 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,0));
449 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
450 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
451 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
452 >      error(CONSISTENCY,"Invalid topology\n");
453 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
454 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
455 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
456 >      error(CONSISTENCY,"Invalid topology\n");
457 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,1));
458 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
459 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
460 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
461 >      error(CONSISTENCY,"Invalid topology\n");
462 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
463 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
464 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
465 >      error(CONSISTENCY,"Invalid topology\n");
466 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,2));
467 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
468 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
469 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
470 >      error(CONSISTENCY,"Invalid topology\n");
471 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
472 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
473 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
474 >      error(CONSISTENCY,"Invalid topology\n");
475 >  }
476 >
477 >  }
478   #endif
479 <            continue;
481 <        }
482 <        t0 = SM_NTH_TRI(sm,id_t0);
483 <        t1 = SM_NTH_TRI(sm,id_t1);
484 <        
485 <        e0 = T_WHICH_V(t0,E_NTH_VERT(e,0));
486 <        T_NTH_NBR(t0,e0) = id_t1;
479 > #endif    
480  
481 <        e1 = T_WHICH_V(t1,E_NTH_VERT(e,1));
482 <        T_NTH_NBR(t1,e1) = id_t0;
483 <    }
484 <    return(test);
481 >    return(TRUE);
482 >
483 > Terror:  
484 >
485 >          error(CONSISTENCY,"smTriangulate():Unable to triangulate\n");
486 >
487   }
488  
489 < eIn_tri(e,t)
490 < int e;
491 < TRI *t;
489 >
490 > void
491 > smTris_swap_edge(sm,t_id,t1_id,e,e1,tn_id,tn1_id)
492 >   SM *sm;
493 >   int t_id,t1_id;
494 >   int e,e1;
495 >   int *tn_id,*tn1_id;
496   {
497 +    S_ID verts[3];
498 +    int enext,eprev,e1next,e1prev;
499 +    TRI *n,*ta,*tb,*t,*t1;
500 +    FVECT p1,p2,p3;
501 +    int ta_id,tb_id;
502 +    /* form new diagonal (e relative to t, and e1 relative to t1)
503 +      defined by quadrilateral formed by t,t1- swap for the opposite diagonal
504 +   */
505 +    t = SM_NTH_TRI(sm,t_id);
506 +    t1 = SM_NTH_TRI(sm,t1_id);
507 +    enext = (e+1)%3;
508 +    eprev = (e+2)%3;
509 +    e1next = (e1+1)%3;
510 +    e1prev = (e1+2)%3;
511 +    verts[e] = T_NTH_V(t,e);
512 +    verts[enext] = T_NTH_V(t,enext);
513 +    verts[eprev] = T_NTH_V(t1,e1);
514 +    ta_id = smAdd_tri(sm,verts[0],verts[1],verts[2],&ta);
515 + #if 0
516 +  fprintf(stderr,"Added tri %d %d %d %d\n",ta_id,T_NTH_V(ta,0),
517 +          T_NTH_V(ta,1), T_NTH_V(ta,2));
518 + #endif
519 +    verts[e1] = T_NTH_V(t1,e1);
520 +    verts[e1next] = T_NTH_V(t1,e1next);
521 +    verts[e1prev] = T_NTH_V(t,e);
522 +    tb_id = smAdd_tri(sm,verts[0],verts[1],verts[2],&tb);
523 + #if 0
524 +  fprintf(stderr,"Added tri %d %d %d %d\n",tb_id,T_NTH_V(tb,0),
525 +          T_NTH_V(tb,1), T_NTH_V(tb,2));
526 + #endif
527 +    /* set the neighbors */
528 +    T_NTH_NBR(ta,e) = T_NTH_NBR(t1,e1next);
529 +    T_NTH_NBR(tb,e1) = T_NTH_NBR(t,enext);
530 +    T_NTH_NBR(ta,enext)= tb_id;
531 +    T_NTH_NBR(tb,e1next)= ta_id;
532 +    T_NTH_NBR(ta,eprev)=T_NTH_NBR(t,eprev);
533 +    T_NTH_NBR(tb,e1prev)=T_NTH_NBR(t1,e1prev);    
534  
535 <  if(T_NTH_V(t,0)==E_NTH_VERT(e,0))
536 <    return(T_NTH_V(t,1)==E_NTH_VERT(e,1)||T_NTH_V(t,2)==E_NTH_VERT(e,1));
537 <  else
538 <    if(T_NTH_V(t,1)==E_NTH_VERT(e,0))
539 <      return(T_NTH_V(t,0)==E_NTH_VERT(e,1)||T_NTH_V(t,2)==E_NTH_VERT(e,1));
540 <    else if(T_NTH_V(t,2)==E_NTH_VERT(e,0))
541 <      return(T_NTH_V(t,0)==E_NTH_VERT(e,1)||T_NTH_V(t,1)==E_NTH_VERT(e,1));
542 <  return(FALSE);
535 >    /* Reset neighbor pointers of original neighbors */
536 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t,enext));
537 >    T_NTH_NBR(n,T_NTH_NBR_PTR(t_id,n)) = tb_id;
538 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t,eprev));
539 >    T_NTH_NBR(n,T_NTH_NBR_PTR(t_id,n)) = ta_id;
540 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,e1next));
541 >    T_NTH_NBR(n,T_NTH_NBR_PTR(t1_id,n)) = ta_id;
542 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,e1prev));
543 >    T_NTH_NBR(n,T_NTH_NBR_PTR(t1_id,n)) = tb_id;
544 >
545 >    smDelete_tri(sm,t_id,t);
546 >    smDelete_tri(sm,t1_id,t1);
547 >
548 > #ifdef DEBUG
549 > #if DEBUG > 1
550 >    if(T_NTH_NBR(ta,0) == T_NTH_NBR(ta,1)  ||
551 >       T_NTH_NBR(ta,1) == T_NTH_NBR(ta,2)  ||
552 >       T_NTH_NBR(ta,0) == T_NTH_NBR(ta,2))
553 >      error(CONSISTENCY,"Invalid topology\n");
554 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,0))) ||
555 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,1))) ||
556 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,2))))
557 >      error(CONSISTENCY,"Invalid topology\n");
558 >    n = SM_NTH_TRI(sm,T_NTH_NBR(ta,0));
559 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
560 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
561 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
562 >      error(CONSISTENCY,"Invalid topology\n");
563 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
564 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
565 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
566 >      error(CONSISTENCY,"Invalid topology\n");
567 >    n = SM_NTH_TRI(sm,T_NTH_NBR(ta,1));
568 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
569 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
570 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
571 >      error(CONSISTENCY,"Invalid topology\n");
572 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
573 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
574 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
575 >      error(CONSISTENCY,"Invalid topology\n");
576 >    n = SM_NTH_TRI(sm,T_NTH_NBR(ta,2));
577 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
578 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
579 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
580 >      error(CONSISTENCY,"Invalid topology\n");
581 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
582 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
583 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
584 >      error(CONSISTENCY,"Invalid topology\n");
585 >    if(T_NTH_NBR(ta,0) == T_NTH_NBR(ta,1)  ||
586 >       T_NTH_NBR(ta,1) == T_NTH_NBR(ta,2)  ||
587 >       T_NTH_NBR(ta,0) == T_NTH_NBR(ta,2))
588 >      error(CONSISTENCY,"Invalid topology\n");
589 >
590 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,0))) ||
591 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,1))) ||
592 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,2))))
593 >      error(CONSISTENCY,"Invalid topology\n");
594 >    n = SM_NTH_TRI(sm,T_NTH_NBR(tb,0));
595 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
596 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
597 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
598 >      error(CONSISTENCY,"Invalid topology\n");
599 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
600 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
601 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
602 >      error(CONSISTENCY,"Invalid topology\n");
603 >    n = SM_NTH_TRI(sm,T_NTH_NBR(tb,1));
604 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
605 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
606 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
607 >      error(CONSISTENCY,"Invalid topology\n");
608 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
609 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
610 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
611 >      error(CONSISTENCY,"Invalid topology\n");
612 >    n = SM_NTH_TRI(sm,T_NTH_NBR(tb,2));
613 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
614 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
615 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
616 >      error(CONSISTENCY,"Invalid topology\n");
617 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
618 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
619 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
620 >      error(CONSISTENCY,"Invalid topology\n");
621 > #endif
622 > #endif
623 >    *tn_id = ta_id;
624 >    *tn1_id = tb_id;
625 >    
626 >    return;
627   }
628 < smFix_edges(sm)
628 >
629 > /* Test the new set of triangles for Delaunay condition. 'Edges' contains
630 >   all of the new edges added. The CCW triangle assoc with each edge is
631 >   tested against the opposite vertex of the CW triangle. If the vertex
632 >   lies inside the circle defined by the CCW triangle- the edge is swapped
633 >   for the opposite diagonal
634 > */
635 > smFixEdges(sm)
636     SM *sm;
637   {
638 <    int e,id_t0,id_t1,e_new,e0,e1,e0_next,e1_next;
639 <    TRI *t0,*t1,*nt0,*nt1;
640 <    int i,id_v0,id_v1,id_v2,id_p,nid_t0,nid_t1;
638 >    int e,t0_id,t1_id,e_new,e0,e1,e0_next,e1_next;
639 >    S_ID v0_id,v1_id,v2_id,p_id;
640 >    int  i,t0_nid,t1_nid;
641      FVECT v0,v1,v2,p,np,v;
642 <    LIST *add,*del;
643 <    
517 <    add = del = NULL;
642 >    TRI *t0,*t1;
643 >
644      FOR_ALL_EDGES(e)
645      {
646 <        id_t0 = E_NTH_TRI(e,0);
647 <        id_t1 = E_NTH_TRI(e,1);
522 <        if((id_t0==SM_INVALID) || (id_t1==SM_INVALID))
523 <        {
646 >        t0_id = E_NTH_TRI(e,0);
647 >        t1_id = E_NTH_TRI(e,1);
648   #ifdef DEBUG
649 <            eputs("smFix_edges: Unassigned edge nbr\n");
649 >        if((t0_id==INVALID) || (t1_id==INVALID))
650 >          error(CONSISTENCY,"smFix_edges: Unassigned edge nbr\n");
651   #endif
652 <            continue;
653 <        }
654 <        t0 = SM_NTH_TRI(sm,id_t0);
655 <        t1 = SM_NTH_TRI(sm,id_t1);
531 <        
532 <        e0 = T_WHICH_V(t0,E_NTH_VERT(e,0));
533 <        e1 = T_WHICH_V(t1,E_NTH_VERT(-e,0));
534 <        e0_next = (e0+2)%3;
535 <        e1_next = (e1+2)%3;
536 <        id_v0 = E_NTH_VERT(e,0);
537 <        id_v1 = E_NTH_VERT(e,1);
538 <        id_v2 = T_NTH_V(t0,e0_next);
539 <        id_p = T_NTH_V(t1,e1_next);
652 >        t0 = SM_NTH_TRI(sm,t0_id);
653 >        t1 = SM_NTH_TRI(sm,t1_id);
654 >        e0 = T_NTH_NBR_PTR(t1_id,t0);
655 >        e1 = T_NTH_NBR_PTR(t0_id,t1);
656  
657 <        smDir(sm,v0,id_v0);
658 <        smDir(sm,v1,id_v1);
659 <        smDir(sm,v2,id_v2);
657 >        v0_id = E_NTH_VERT(e,0);
658 >        v1_id = E_NTH_VERT(e,1);
659 >        v2_id = T_NTH_V(t0,e0);
660 >        p_id = T_NTH_V(t1,e1);
661 >
662 >        smDir(sm,v0,v0_id);
663 >        smDir(sm,v1,v1_id);
664 >        smDir(sm,v2,v2_id);
665          
666 <        VCOPY(p,SM_NTH_WV(sm,id_p));    
667 <        VSUB(p,p,SM_VIEW_CENTER(sm));
668 <        if(point_in_cone(p,v0,v1,v2))
666 >        VSUB(p,SM_NTH_WV(sm,p_id),SM_VIEW_CENTER(sm));
667 >        normalize(p);
668 >        if(pt_in_cone(p,v2,v1,v0))
669          {
670 <            smTris_swap_edge(sm,id_t0,id_t1,e0,e1,&nid_t0,&nid_t1,&add,&del);
671 <            
672 <            nt0 = SM_NTH_TRI(sm,nid_t0);
673 <            nt1 = SM_NTH_TRI(sm,nid_t1);
670 >           smTris_swap_edge(sm,t0_id,t1_id,e0,e1,&t0_nid,&t1_nid);
671 >           /* Adjust the triangle pointers of the remaining edges to be
672 >              processed
673 >            */
674              FOR_ALL_EDGES_FROM(e,i)
675              {
676 <              if(E_NTH_TRI(i,0)==id_t0 || E_NTH_TRI(i,0)==id_t1)
676 >              if(E_NTH_TRI(i,0)==t0_id || E_NTH_TRI(i,0)==t1_id)
677                {
678 <                if(eIn_tri(i,nt0))
679 <                  SET_E_NTH_TRI(i,0,nid_t0);
678 >                if(eIn_tri(i,SM_NTH_TRI(sm,t0_nid)))
679 >                  SET_E_NTH_TRI(i,0,t0_nid);
680                  else
681 <                  SET_E_NTH_TRI(i,0,nid_t1);
681 >                  SET_E_NTH_TRI(i,0,t1_nid);
682                }
683  
684 <              if(E_NTH_TRI(i,1)==id_t0 || E_NTH_TRI(i,1)==id_t1)
684 >              if(E_NTH_TRI(i,1)==t0_id || E_NTH_TRI(i,1)==t1_id)
685                {
686 <                if(eIn_tri(i,nt0))
687 <                  SET_E_NTH_TRI(i,1,nid_t0);
686 >                if(eIn_tri(i,SM_NTH_TRI(sm,t0_nid)))
687 >                  SET_E_NTH_TRI(i,1,t0_nid);
688                  else
689 <                  SET_E_NTH_TRI(i,1,nid_t1);
689 >                  SET_E_NTH_TRI(i,1,t1_nid);
690                }
691 + #ifdef DEBUG
692 +              if(E_NTH_TRI(i,1) == E_NTH_TRI(i,0) )
693 +                error(CONSISTENCY,"invalid edge\n");
694 + #endif
695              }
696 <            id_t0 = nid_t0;
697 <            id_t1 = nid_t1;
696 >            t0_id = t0_nid;
697 >            t1_id = t1_nid;
698              e_new = eNew_edge();
699 <            SET_E_NTH_VERT(e_new,0,id_p);
700 <            SET_E_NTH_VERT(e_new,1,id_v2);
701 <            SET_E_NTH_TRI(e_new,0,id_t0);
702 <            SET_E_NTH_TRI(e_new,1,id_t1);
699 >            SET_E_NTH_VERT(e_new,0,p_id);
700 >            SET_E_NTH_VERT(e_new,1,v2_id);
701 >            SET_E_NTH_TRI(e_new,0,t0_id);
702 >            SET_E_NTH_TRI(e_new,1,t1_id);
703 > #ifdef DEBUG
704 >              if(E_NTH_TRI(i,1) == E_NTH_TRI(i,0) )
705 >                error(CONSISTENCY,"invalid edge\n");
706 > #endif
707          }
708      }
580 –    smUpdate_locator(sm,add,del);
709   }
710  
711 < int
712 < smMesh_remove_vertex(sm,id)
713 <   SM *sm;
714 <   int id;
711 >
712 > smDelete_samp(sm,s_id)
713 > SM *sm;
714 > S_ID s_id;
715   {
716 <    int tri;
717 <    LIST *elist;
590 <    int cnt,debug;
591 <    /* generate list of vertices that form the boundary of the
592 <       star polygon formed by vertex id and all of its adjacent
593 <       triangles
594 <     */
595 <    eClear_edges();
596 <    elist = smVertex_star_polygon(sm,id);
597 <    if(!elist)
598 <       return(FALSE);
716 >  QUADTREE qt;
717 >  S_ID *os;
718  
719 <    /* Triangulate spherical polygon */
720 <    smTriangulate(sm,elist);
719 >  /* Mark as free */
720 >  smUnalloc_samp(sm,s_id);
721  
722 <
723 <    /* Fix up new triangles to be Delaunay */
724 <    smFix_edges(sm);
725 <
726 <    return(TRUE);
722 > #ifdef DEBUG  
723 >  SM_NTH_VERT(sm,s_id) = INVALID;
724 >  /* fprintf(stderr,"deleting sample %d\n",s_id); */
725 > #endif
726 >  /* remove from its set */
727 >  qt = SM_S_NTH_QT(sm,s_id);
728 >  os = qtqueryset(qt);
729 >  deletuelem(os, s_id); /* delete obj from unsorted os, no questions */
730   }
731 <  
732 < /* Remove point from samples, and from mesh. Delete any triangles
733 <   adjacent to the point and re-triangulate the hole
612 <   Return TRUE is point found , FALSE otherwise
613 < */
731 > /* Remove vertex "id" from the mesh- and retriangulate the resulting
732 >   hole: Returns TRUE if successful, FALSE otherwise.
733 > */
734   int
735 < smDelete_point(sm,id)
736 < SM *sm;
737 < int id;
735 > smRemoveVertex(sm,id)
736 >   SM *sm;
737 >   S_ID id;
738   {
739 +    LIST *b_list;
740 +    /* generate list of edges that form the boundary of the
741 +       polygon formed by the triangles adjacent to vertex 'id'*/
742 +    b_list = smVertexStar(sm,id);
743 + #if 0
744 + {int i;
745 +    eputs("\n\n");
746 +    for(i=1;i<=Ecnt;i++)
747 +      fprintf(stderr,"%d verts %d %d tris  %d %d\n",
748 +              i,Edges[i].verts[0],Edges[i].verts[1],
749 +              Edges[i].tris[0],Edges[i].tris[1]);
750 + }
751 + #endif
752  
753 <    /* Remove the corresponding vertex from the mesh */
754 <    smMesh_remove_vertex(sm,id);
755 <    /* Free the sample point */
756 <    smDelete_sample(sm,id);
753 >    /* Triangulate polygonal hole  */
754 >    smTriangulate(sm,id,b_list);
755 >    
756 >    /* Fix up new triangles to be Delaunay*/
757 >
758 >    smFixEdges(sm);
759 >    smDelete_samp(sm,id);
760 >    eClear_edges();
761      return(TRUE);
762   }
763 +  
764  
765  
766 <
766 >
767  
768  
769  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines