--- ray/src/hd/sm_del.c	1998/09/11 11:52:25	3.3
+++ ray/src/hd/sm_del.c	1999/06/10 15:22:22	3.12
@@ -8,487 +8,481 @@ static char SCCSid[] = "$SunId$ SGI";
  *  sm_del.c
  */
 #include "standard.h"
-
+#include "sm_flag.h"
 #include "sm_list.h"
 #include "sm_geom.h"
 #include "sm.h"
 
-static EDGE Edges[MAX_EDGES];
+#define MAX_EDGE_INIT 100
+static int Max_edges= MAX_EDGE_INIT;
+static EDGE *Edges=NULL;
 static int Ecnt=0;
 
-int
-remove_tri(qtptr,fptr,t_id)
-   QUADTREE *qtptr;
-   int *fptr;
-   int t_id;
-{
-    OBJECT tset[QT_MAXSET+1];	
-    int n;
-
-   if(QT_IS_EMPTY(*qtptr))
-     return(FALSE);
-   /* remove id from set */
-      else
-      {
-	 if(!qtinset(*qtptr,t_id))
-	   return(FALSE);
-	 n = QT_SET_CNT(qtqueryset(*qtptr))-1;
-	 *qtptr = qtdelelem(*qtptr,t_id);
-	  if(n  == 0)
-	     (*fptr) |= QT_COMPRESS;
-	if(!QT_FLAG_FILL_TRI(*fptr))
-	  (*fptr)++;
-      }
-    return(TRUE);
-}
-
-int
-remove_tri_compress(qtptr,q0,q1,q2,t0,t1,t2,n,arg,t_id)
-QUADTREE *qtptr;
-FVECT q0,q1,q2;
-FVECT t0,t1,t2;
-int n;
-int *arg;
-int t_id;
-{
-    int f = 0;
-    /* NOTE compress */
-    return(remove_tri(qtptr,&f,t_id));
-}
-
-
-
-int
-stDelete_tri(st,t_id,t0,t1,t2)
-   STREE *st;
-   int t_id;
-   FVECT t0,t1,t2;
-{
-    int f;
-    FVECT dir;
-    
-  /* First add all of the leaf cells lying on the triangle perimeter:
-     mark all cells seen on the way
-   */
-    ST_CLEAR_FLAGS(st);
-    f = 0;
-    VSUB(dir,t1,t0);
-    stTrace_edge(st,t0,dir,1.0,remove_tri,&f,t_id);
-    VSUB(dir,t2,t1);
-    stTrace_edge(st,t1,dir,1.0,remove_tri,&f,t_id);
-    VSUB(dir,t0,t2);
-    stTrace_edge(st,t2,dir,1.0,remove_tri,&f,t_id);
-    /* Now visit interior */
-    if(QT_FLAG_FILL_TRI(f) || QT_FLAG_UPDATE(f))
-       stVisit_tri_interior(st,t0,t1,t2,remove_tri_compress,&f,t_id);
-}
-
-
-smLocator_remove_tri(sm,t_id)
+smFree_tri(sm,id,t)
 SM *sm;
-int t_id;
-{
-  STREE *st;
-  char found;
-  TRI *t;
-  FVECT v0,v1,v2;
-
-  st = SM_LOCATOR(sm);
-
-  t = SM_NTH_TRI(sm,t_id);
-
-  VSUB(v0,SM_T_NTH_WV(sm,t,0),SM_VIEW_CENTER(sm));
-  VSUB(v1,SM_T_NTH_WV(sm,t,1),SM_VIEW_CENTER(sm));
-  VSUB(v2,SM_T_NTH_WV(sm,t,2),SM_VIEW_CENTER(sm));
-  found = stUpdate_tri(st,t_id,v0,v1,v2,remove_tri,remove_tri_compress);
-  return(found);
-}
-
-smFree_tri(sm,id)
-SM *sm;
 int id;
+TRI *t;
 {
-  TRI *tri;
-
-  tri = SM_NTH_TRI(sm,id);
   /* Add to the free_list */
-  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
      has been added- before the deletion: replacing
      the old tri- and therefore dont need to dereference any pointers
      to id because the vertices can no longer
      point to tri id as being the first triangle pointer
   */
-  if(!SM_IS_NTH_T_BASE(sm,t_id))
-  {
-    SM_NUM_TRIS(sm)--;
-    if(SM_IS_NTH_T_NEW(sm,t_id))
-      smNew_tri_cnt--;
-  }
-  smClear_tri_flags(sm,t_id);
+  SM_CLR_NTH_T_ACTIVE(sm,t_id);
+  smFree_tri(sm,t_id,t);
 
-  smFree_tri(sm,t_id);
-
 }
 
+int
+eNew_edge()
+{
+  if(!Edges)
+    if(!(Edges = (EDGE *)realloc(NULL,(Max_edges+1)*sizeof(EDGE))))
+      goto memerr;
 
+  if(Ecnt >= Max_edges)
+    {
+#ifdef DEBUG
+      if(Max_edges > 10000)
+	error(CONSISTENCY,"Too many edges in vertex loop\n");
+#endif
+      Max_edges += 100;
+      if(!(Edges = (EDGE *)realloc(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);
 
-LIST
-*smVertex_star_polygon(sm,id)
+memerr:
+  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)
 SM *sm;
 int id;
 {
     TRI *tri,*t_next;
-    LIST *elist,*end,*tlist;
-    int t_id,v_next,t_next_id;
-    int e;
+    LIST *elist,*end;
+    int e,t_id,v_next,t_next_id,b_id,v_id,t_last_id;
 
     elist = end =  NULL;
+
     /* Get the first triangle adjacent to vertex id */
     t_id = SM_NTH_VERT(sm,id);
     tri = SM_NTH_TRI(sm,t_id);
 
-    
-    if((e = eNew_edge()) == SM_INVALID)
-    {
-#ifdef DEBUG
-	eputs("smVertex_star_polygon():Too many edges\n");
-#endif
-	return(NULL);
-    }
+    e = eNew_edge();
+    /* Get the  next vertex on the polygon boundary */
+    v_id = T_WHICH_V(tri,id);
+    b_id = (v_id + 1)%3;
+    /* Create an edge */
+    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);
+    SET_E_NTH_VERT(e,1,v_next);
+  
     elist = add_data_to_circular_list(elist,&end,e);
-    v_next = (T_WHICH_V(tri,id)+1)%3;
-    SET_E_NTH_VERT(e,0,T_NTH_V(tri,v_next));
-    SET_E_NTH_TRI(e,0,SM_INVALID);
-    SET_E_NTH_TRI(e,1,T_NTH_NBR(tri,v_next));
-    v_next = (T_WHICH_V(tri,id)+2)%3;
-    SET_E_NTH_VERT(e,1,T_NTH_V(tri,v_next));
-
-    
     t_next_id = t_id;
     t_next = tri;
+    t_last_id = -1;
+ 
+    /* Create a set to hold all of the triangles for deletion later */
 
-    tlist = push_data(NULL,t_id);
+    while((t_next_id = T_NTH_NBR(t_next,b_id)) != t_id)
+    {
+      e = eNew_edge();
+      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);
+      b_id = (v_id + 1)%3;
+      SET_E_NTH_TRI(e,1,T_NTH_NBR(t_next,v_id));
+      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);
 
-    while((t_next_id = smTri_next_ccw_nbr(sm,t_next,id)) != t_id)
-    {	
-	if((e = eNew_edge()) == SM_INVALID)
-	{
-#ifdef DEBUG
-	    eputs("smVertex_star_polygon():Too many edges\n");
-#endif
-	    return(NULL);
-	}
-	elist = add_data_to_circular_list(elist,&end,e);
-	t_next = SM_NTH_TRI(sm,t_next_id);
-	v_next = (T_WHICH_V(t_next,id)+1)%3;
-	SET_E_NTH_VERT(e,0,T_NTH_V(t_next,v_next));
-	SET_E_NTH_TRI(e,0,SM_INVALID);
-	SET_E_NTH_TRI(e,1,T_NTH_NBR(t_next,v_next));
-	v_next = (T_WHICH_V(t_next,id)+2)%3;
-	SET_E_NTH_VERT(e,1,T_NTH_V(t_next,v_next));
-	tlist = push_data(tlist,t_next_id);
     }
-    while(tlist)
-    {
-	t_id = (int)pop_list(&tlist);
-	/* first remove from point location structure */
-	smLocator_remove_tri(sm,t_id);
-	smDelete_tri(sm,t_id);
-    }
-    return(elist);
+    smDelete_tri(sm,t_last_id,t_next);
+    smDelete_tri(sm,t_id,tri); 
+   return(elist);
 }
 
 int
-smEdge_intersect_polygon(sm,v0,v1,l)
+smTriangulate_add_tri(sm,id0,id1,id2,e0,e1,e2ptr)
 SM *sm;
-FVECT v0,v1;
-LIST *l;
+int id0,id1,id2,e0,e1,*e2ptr;
 {
-    FVECT e0,e1;
-    int e,id_e0,id_e1;
-    LIST *el,*eptr;
-    
-    /* Test the edges in l against v0v1 to see if v0v1 intersects
-       any other edges
-     */
-    
-    el = l;
+  int t_id,e2;
+  TRI *t;
 
-    while(el)
-    {
-      e = (int)LIST_DATA(el);
-      id_e0 = E_NTH_VERT(e,0);
-      id_e1 = E_NTH_VERT(e,1);
-      /* NOTE: DO these need to be normalized? Just subtract center? */
-      smDir(sm,e0,id_e0);
-      smDir(sm,e1,id_e1);
-      if(sedge_intersect(v0,v1,e0,e1))
-	return(TRUE);
+  t_id = smAdd_tri(sm,id0,id1,id2,&t);
+  if(*e2ptr == 0)
+  {
+    e2 = eNew_edge();
+    SET_E_NTH_VERT(e2,0,id2);
+    SET_E_NTH_VERT(e2,1,id0);
+  }
+  else
+    e2 = *e2ptr;
+  /* set appropriate tri for each edge*/ 
+  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);
 
-      el = LIST_NEXT(el);
-      if(el == l)
-	break;
-    }
-    return(FALSE);
 }
 
-int
-smFind_next_convex_vertex(sm,id0,id1,v0,v1,l)
-   SM *sm;
-   int id0,id1;
-   FVECT v0,v1;
-   LIST *l;
+int 
+smTriangulate_quad(sm,l)
+     SM *sm;
+     LIST *l;
 {
-    int e,id;
-    LIST *el;
-    FVECT v;
+  int e1,e2,e3,e4,e,t_id,id0,id1,id2,id3;
+  FVECT v0,v1,v2,v3,n;
+  double dp;
+  TRI *tri;
+  LIST *lptr;
 
-    /* starting with the end of edge at head of l, search sequentially for
-      vertex v such that v0v1v is a convex angle, and the edge v1v does
-      not intersect any other edges
-   */
-    id = SM_INVALID;
-    el = l;
-    while(id != id0)
-    {
-	e = (int)LIST_DATA(el);
-	id = E_NTH_VERT(e,1);
-
-	smDir(sm,v,id); 
-
-	if(convex_angle(v0,v1,v) && !smEdge_intersect_polygon(sm,v1,v,l))
-	   return(id);
-	      
-	el = LIST_NEXT(el);
-	if(el == l)
-	   break;
-    }
-    return(SM_INVALID);
-}
-
-int
-split_edge_list(id0,id_new,l,lnew)
-int id0,id_new;
-LIST **l,**lnew;
-{
-    LIST *list,*lptr,*end;
-    int e,e1,e2,new_e;
-
-    e2 = SM_INVALID;
-    list = lptr = *l;
-
-    if((new_e = eNew_edge())==SM_INVALID)
-     {
 #ifdef DEBUG
-	    eputs("split_edge_list():Too many edges\n");
+  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
-	 return(FALSE);
-     }
-    SET_E_NTH_VERT(new_e,0,id0);
-    SET_E_NTH_VERT(new_e,1,id_new);
-    SET_E_NTH_TRI(new_e,0,SM_INVALID);
-    SET_E_NTH_TRI(new_e,1,SM_INVALID);
-    
-    while(e2 != id_new)
-    {
-	lptr = LIST_NEXT(lptr);
-	e = (int)LIST_DATA(lptr);
-	e2 = E_NTH_VERT(e,1);
-	if(lptr == list)
-	{
-#ifdef DEBUG	    
-	  eputs("split_edge_list():cant find vertex\n");
-#endif
-	  *lnew = NULL;
-	  return(FALSE);
-	}
+  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);
 
+  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);
+      }
     }
-    end = lptr;
-    lptr =  LIST_NEXT(lptr);
-    list = add_data_to_circular_list(list,&end,-new_e);
-    *lnew = list;
 
-    /* now follow other cycle */
-
-    list = lptr;
-    e2 = SM_INVALID;
-    while(e2 != id0)
-    {
-	lptr = LIST_NEXT(lptr);
-	e = (int)LIST_DATA(lptr);
-	e2 = E_NTH_VERT(e,1);
-	if(lptr == list)
-	{
-#ifdef DEBUG	    
-	  eputs("split_edge_list():cant find intial vertex\n");
+  }
+#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
-	  *l = NULL;
-	  return(FALSE);
-	}
-
-    }
-    end = lptr;
-    list = add_data_to_circular_list(list,&end,new_e);
-    *l = list;
-    return(TRUE);
+#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);
 }
 
+/* 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,
+   assumes polygon star relative to id
+*/
 
 int
-smTriangulate_convex(sm,plist)
+smTriangulate(sm,id,plist)
 SM *sm;
+int id;
 LIST *plist;
 {
-    TRI *tri;
-    int t_id,e_id0,e_id1,e_id2;
-    int v_id0,v_id1,v_id2;
-    LIST *lptr;
-    int cnt;
+    LIST *l,*prev,*t;
+    FVECT v0,v1,v2,n,p;
+    int is_tri,cut,t_id,id0,id1,id2,e2,e1,enew;
+    double dp1,dp2;
 
-    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);
-	/* form a triangle for each triple of with v0 as base of star */
-	t_id = smAdd_tri(sm,v_id0,v_id1,v_id2,&tri);
-	smLocator_add_tri(sm,t_id,v_id0,v_id1,v_id2);
-	/* add which pointer?*/
-
-	lptr = LIST_NEXT(lptr);
-
-	if(LIST_NEXT(lptr) != plist)
-	{
-	    e_id2 = eNew_edge();
-	    SET_E_NTH_VERT(e_id2,0,v_id2);
-	    SET_E_NTH_VERT(e_id2,1,v_id0);
-	}
-	else
-	   e_id2 = (int)LIST_DATA(lptr);
-	 
-	/* set appropriate tri for each edge*/ 
-	SET_E_NTH_TRI(e_id0,0,t_id);
-	SET_E_NTH_TRI(e_id1,0,t_id);
-	SET_E_NTH_TRI(e_id2,0,t_id);
-
-	e_id0 = -e_id2;
-    }
-    
-    free_list(plist);
-    return(TRUE);
-}
-int
-smTriangulate_elist(sm,plist)
-SM *sm;
-LIST *plist;
-{
-    LIST *l,*el1;
-    FVECT v0,v1,v2;
-    int id0,id1,id2,e,id_next;
-    char flipped;
-    int done;
-
+    VSUB(p,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
+    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));   
+    normalize(v1);
     while(l)
     {
-	/* get v0,v1,v2 */
-      e = (int)LIST_DATA(l);
-      id0 = E_NTH_VERT(e,0);
-      id1 = E_NTH_VERT(e,1);
       l = LIST_NEXT(l);
-      e = (int)LIST_DATA(l);
-      id2 = E_NTH_VERT(e,1);
-
-      smDir(sm,v0,id0);
-      smDir(sm,v1,id1);
-      smDir(sm,v2,id2);
-      /* determine if convex (left turn), or concave(right turn) angle */
-      if(convex_angle(v0,v1,v2))
+      /* Get v2 */
+      e2 = (int)LIST_DATA(l);
+      id2 = E_NTH_VERT(e2,1);
+      VSUB(v2,SM_NTH_WV(sm,id2),SM_VIEW_CENTER(sm));
+      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
+       */
+      VCROSS(n,v0,v2);
+      normalize(n);
+      dp1 = DOT(n,p);      
+      if(dp1  <= 0.0)
       {
-	if(l == plist)
-	  break;
-	else
-	  continue;
+	  /* 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
+        */
+	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;
+	    /* 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;
+	}
       }
-      /* if concave: add edge and recurse on two sub polygons */
-      id_next = smFind_next_convex_vertex(sm,id0,id1,v0,v1,LIST_NEXT(l));
-      if(id_next == SM_INVALID)
+      VCOPY(v0,v1);
+      VCOPY(v1,v2);
+      id0 = id1;
+      id1 = id2;
+      e1 = e2;	
+      /* check if gone around circular list without adding any
+	 triangles: prevent infinite loop */
+      if(l == plist)
       {
-#ifdef DEBUG
-	  eputs("smTriangulate_elist():Unable to find convex vertex\n");
-#endif
-	  return(FALSE);
+	if(LIST_NEXT(LIST_NEXT(l)) == prev)
+	  {/* Check if have a triangle */
+	    is_tri = TRUE;	
+	    break;
+	  }
+	if(!cut)
+	  break;
+        cut = FALSE;
       }
-      /* add edge */
-      el1 = NULL;
-      /* Split edge list l into two lists: one from id1-id_next-id1,
-	 and the next from id2-id_next-id2
-      */
-      split_edge_list(id1,id_next,&l,&el1);
-      /* Recurse and triangulate the two edge lists */
-      done = smTriangulate_elist(sm,l);
-      if(done)
-	done = smTriangulate_elist(sm,el1);
-      return(done);
+      prev = l;
     }
-    done = smTriangulate_convex(sm,plist);
-    return(done);
-}
-
-int
-smTriangulate(sm,plist)
-SM *sm;
-LIST *plist;
-{
-    int e,id_t0,id_t1,e0,e1;
-    TRI *t0,*t1;
-    int test;
-    
-    test = smTriangulate_elist(sm,plist);
-
-    if(!test)
-       return(test);
-    FOR_ALL_EDGES(e)
+    if(is_tri)
     {
-	id_t0 = E_NTH_TRI(e,0);
-	id_t1 = E_NTH_TRI(e,1);
-	if((id_t0==SM_INVALID) || (id_t1==SM_INVALID))
-	{
-#ifdef DEBUG
-	   eputs("smTriangulate(): Unassigned edge neighbor\n");
+      l = LIST_NEXT(l);
+      enew = (int)LIST_DATA(l);
+      t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew);
+      free_list(l);
+    }
+    else
+      if(!smTriangulate_quad(sm,l))
+	goto Terror;
+#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
-	    continue;
-	}
-	t0 = SM_NTH_TRI(sm,id_t0);
-	t1 = SM_NTH_TRI(sm,id_t1);
+
+    /* Set triangle adjacencies based on edge adjacencies */ 
+    FOR_ALL_EDGES(enew)
+    {
+      id0 = E_NTH_TRI(enew,0);
+      id1 = E_NTH_TRI(enew,1);
 	
-	e0 = T_WHICH_V(t0,E_NTH_VERT(e,0));
-	T_NTH_NBR(t0,e0) = id_t1;
+      e1 = (T_WHICH_V(SM_NTH_TRI(sm,id0),E_NTH_VERT(enew,0))+2)%3;
+      T_NTH_NBR(SM_NTH_TRI(sm,id0),e1) = id1;
+      
+      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;
+    }
+#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
 
-	e1 = T_WHICH_V(t1,E_NTH_VERT(e,1));
-	T_NTH_NBR(t1,e1) = id_t0;
+#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");
     }
-    return(test);
+    }
+#endif
+#endif
+    return(TRUE);
+
+Terror:	  
+#ifdef DEBUG
+	  error(CONSISTENCY,"smTriangulate():Unable to triangulate\n");
+#endif
 }
 
 eIn_tri(e,t)
@@ -503,129 +497,286 @@ TRI *t;
       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);
 }
-smFix_edges(sm)
+
+
+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;
 {
-    int e,id_t0,id_t1,e_new,e0,e1,e0_next,e1_next;
-    TRI *t0,*t1,*nt0,*nt1;
-    int i,id_v0,id_v1,id_v2,id_p,nid_t0,nid_t1;
-    FVECT v0,v1,v2,p,np,v;
-    LIST *add,*del;
+    int verts[3],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);    
+
+    /* 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;
+
+    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;
     
-    add = del = NULL;
+    return;
+}
+
+/* Test the new set of triangles for Delaunay condition. 'Edges' contains
+   all of the new edges added. The CCW triangle assoc with each edge is
+   tested against the opposite vertex of the CW triangle. If the vertex
+   lies inside the circle defined by the CCW triangle- the edge is swapped 
+   for the opposite diagonal
+*/
+smFixEdges(sm)
+   SM *sm;
+{
+    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;
+    FVECT v0,v1,v2,p,np,v;
+    TRI *t0,*t1;
+
     FOR_ALL_EDGES(e)
     {
-	id_t0 = E_NTH_TRI(e,0);
-	id_t1 = E_NTH_TRI(e,1);
-	if((id_t0==SM_INVALID) || (id_t1==SM_INVALID))
-	{
+	t0_id = E_NTH_TRI(e,0);
+	t1_id = E_NTH_TRI(e,1);
 #ifdef DEBUG
-	    eputs("smFix_edges: Unassigned edge nbr\n");
+	if((t0_id==INVALID) || (t1_id==INVALID))
+	  error(CONSISTENCY,"smFix_edges: Unassigned edge nbr\n");
 #endif
-	    continue;
-	}
-	t0 = SM_NTH_TRI(sm,id_t0);
-	t1 = SM_NTH_TRI(sm,id_t1);
-	
-	e0 = T_WHICH_V(t0,E_NTH_VERT(e,0));
-	e1 = T_WHICH_V(t1,E_NTH_VERT(-e,0));
-	e0_next = (e0+2)%3;
-	e1_next = (e1+2)%3;
-	id_v0 = E_NTH_VERT(e,0);
-	id_v1 = E_NTH_VERT(e,1);
-	id_v2 = T_NTH_V(t0,e0_next);
-	id_p = T_NTH_V(t1,e1_next);
+	t0 = SM_NTH_TRI(sm,t0_id);
+	t1 = SM_NTH_TRI(sm,t1_id);
+	e0 = T_NTH_NBR_PTR(t1_id,t0);
+	e1 = T_NTH_NBR_PTR(t0_id,t1);
 
-	smDir(sm,v0,id_v0);
-	smDir(sm,v1,id_v1);
-	smDir(sm,v2,id_v2);
+	v0_id = E_NTH_VERT(e,0);
+	v1_id = E_NTH_VERT(e,1);
+	v2_id = T_NTH_V(t0,e0);
+	p_id = T_NTH_V(t1,e1);
+
+	smDir(sm,v0,v0_id);
+	smDir(sm,v1,v1_id);
+	smDir(sm,v2,v2_id);
 	
-	VCOPY(p,SM_NTH_WV(sm,id_p));	
-	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,id_t0,id_t1,e0,e1,&nid_t0,&nid_t1,&add,&del);
-	    
-	    nt0 = SM_NTH_TRI(sm,nid_t0);
-	    nt1 = SM_NTH_TRI(sm,nid_t1);
+	   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
+            */
 	    FOR_ALL_EDGES_FROM(e,i)
 	    {
-	      if(E_NTH_TRI(i,0)==id_t0 || E_NTH_TRI(i,0)==id_t1)
+	      if(E_NTH_TRI(i,0)==t0_id || E_NTH_TRI(i,0)==t1_id)
 	      {
-		if(eIn_tri(i,nt0))
-		  SET_E_NTH_TRI(i,0,nid_t0);
+		if(eIn_tri(i,SM_NTH_TRI(sm,t0_nid)))
+		  SET_E_NTH_TRI(i,0,t0_nid);
 		else
-		  SET_E_NTH_TRI(i,0,nid_t1);
+		  SET_E_NTH_TRI(i,0,t1_nid);
 	      }
 
-	      if(E_NTH_TRI(i,1)==id_t0 || E_NTH_TRI(i,1)==id_t1)
+	      if(E_NTH_TRI(i,1)==t0_id || E_NTH_TRI(i,1)==t1_id)
 	      {
-		if(eIn_tri(i,nt0))
-		  SET_E_NTH_TRI(i,1,nid_t0);
+		if(eIn_tri(i,SM_NTH_TRI(sm,t0_nid)))
+		  SET_E_NTH_TRI(i,1,t0_nid);
 		else
-		  SET_E_NTH_TRI(i,1,nid_t1);
+		  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
 	    }
-	    id_t0 = nid_t0;
-	    id_t1 = nid_t1;
+	    t0_id = t0_nid;
+	    t1_id = t1_nid;
 	    e_new = eNew_edge();
-	    SET_E_NTH_VERT(e_new,0,id_p);
-	    SET_E_NTH_VERT(e_new,1,id_v2);
-	    SET_E_NTH_TRI(e_new,0,id_t0);
-	    SET_E_NTH_TRI(e_new,1,id_t1);
+	    SET_E_NTH_VERT(e_new,0,p_id);
+	    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
 	}
     }
-    smUpdate_locator(sm,add,del);
 }
 
+
+smDelete_samp(sm,s_id)
+SM *sm;
+int s_id;
+{
+  QUADTREE qt;
+  OBJECT *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
-smMesh_remove_vertex(sm,id)
+smRemoveVertex(sm,id)
    SM *sm;
    int id;
 {
-    int tri;
-    LIST *elist;
-    int cnt,debug;
-    /* generate list of vertices that form the boundary of the
-       star polygon formed by vertex id and all of its adjacent
-       triangles
-     */
-    eClear_edges();
-    elist = smVertex_star_polygon(sm,id);
-    if(!elist)
-       return(FALSE);
+    LIST *b_list;
+    /* generate list of edges that form the boundary of the
+       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 spherical polygon */
-    smTriangulate(sm,elist);
+    /* Triangulate polygonal hole  */
+    smTriangulate(sm,id,b_list);
+    
+    /* Fix up new triangles to be Delaunay*/
 
-
-    /* Fix up new triangles to be Delaunay */
-    smFix_edges(sm);
-
+    smFixEdges(sm);
+    smDelete_samp(sm,id);
+    eClear_edges();
     return(TRUE);
 }
    
-/* Remove point from samples, and from mesh. Delete any triangles
-   adjacent to the point and re-triangulate the hole
-   Return TRUE is point found , FALSE otherwise
-*/
-int
-smDelete_point(sm,id)
-SM *sm;
-int id;
-{
 
-    /* Remove the corresponding vertex from the mesh */
-    smMesh_remove_vertex(sm,id);
-    /* Free the sample point */
-    smDelete_sample(sm,id);
-    return(TRUE);
-}
 
-
-
+