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.12 by gwlarson, Thu Jun 10 15:22:22 1999 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines