--- ray/src/hd/sm_del.c	1998/10/06 18:16:53	3.6
+++ ray/src/hd/sm_del.c	1999/06/10 15:22:22	3.12
@@ -11,75 +11,33 @@ static char SCCSid[] = "$SunId$ SGI";
 #include "sm_flag.h"
 #include "sm_list.h"
 #include "sm_geom.h"
-#include "sm_qtree.h"
-#include "sm_stree.h"
 #include "sm.h"
 
-static int Max_edges=200;
+#define MAX_EDGE_INIT 100
+static int Max_edges= MAX_EDGE_INIT;
 static EDGE *Edges=NULL;
 static int Ecnt=0;
 
-#define remove_tri_compress remove_tri
-remove_tri(qtptr,fptr,tptr)
-   QUADTREE *qtptr;
-   int *fptr;
-   int *tptr;
-{
-    int n;
-
-    if(QT_IS_EMPTY(*qtptr))
-      return;
-    if(QT_LEAF_IS_FLAG(*qtptr))
-      return;
-
-    n = QT_SET_CNT(qtqueryset(*qtptr))-1;
-    *qtptr = qtdelelem(*qtptr,*tptr);
-    if(n  == 0)
-      (*fptr) |= QT_COMPRESS;
-    if(!QT_FLAG_FILL_TRI(*fptr))
-      (*fptr)++;
-}
-
-
-smLocator_remove_tri(sm,t_id,v0_id,v1_id,v2_id)
+smFree_tri(sm,id,t)
 SM *sm;
-int t_id;
-int v0_id,v1_id,v2_id;
-{
-  STREE *st;
-  FVECT v0,v1,v2;
-  
-  st = SM_LOCATOR(sm);
-
-  VSUB(v0,SM_NTH_WV(sm,v0_id),SM_VIEW_CENTER(sm));
-  VSUB(v1,SM_NTH_WV(sm,v1_id),SM_VIEW_CENTER(sm));
-  VSUB(v2,SM_NTH_WV(sm,v2_id),SM_VIEW_CENTER(sm));
-
-  qtClearAllFlags();
-  
-  stApply_to_tri(st,v0,v1,v2,remove_tri,remove_tri_compress,&t_id);
-
-}
-
-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
@@ -88,18 +46,11 @@ int t_id;
      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_SAMPLE_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()
 {
@@ -109,333 +60,93 @@ eNew_edge()
 
   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);
 
 memerr:
-  error(SYSTEM,"eNew_edge(): Unable to allocate memory");
+  error(SYSTEM,"eNew_edge): Unable to allocate memory");
 }
 
+/* Return list of edges defining polygon formed by boundary of triangles 
+adjacent to id. Return set of triangles adjacent to id to delete in delptr
+*/
 LIST 
-*smVertex_star_polygon(sm,id,delptr)
+*smVertexStar(sm,id)
 SM *sm;
 int id;
-QUADTREE *delptr;
 {
     TRI *tri,*t_next;
     LIST *elist,*end;
-    int t_id,v_next,t_next_id;
-    int e;
-    OBJECT del_set[2];
+    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()) == INVALID)
-      return(NULL);
-
-    v_next = (T_WHICH_V(tri,id)+1)%3;
-    SET_E_NTH_VERT(e,0,T_NTH_V(tri,v_next));
+    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_next));
-    v_next = (T_WHICH_V(tri,id)+2)%3;
-    SET_E_NTH_VERT(e,1,T_NTH_V(tri,v_next));
+    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);
-
     t_next_id = t_id;
     t_next = tri;
+    t_last_id = -1;
+ 
+    /* Create a set to hold all of the triangles for deletion later */
 
-    del_set[0] =1; del_set[1] = t_id;
-    *delptr = qtnewleaf(del_set);
-
-    while((t_next_id = T_NTH_NBR(t_next,v_next)) != t_id)
-    {	
-      if((e = eNew_edge()) == INVALID)
-	return(NULL);
-
+    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);
-      v_next = (T_WHICH_V(t_next,id)+1)%3;
-
-      SET_E_NTH_VERT(e,0,T_NTH_V(t_next,v_next));
+      t_last_id = t_next_id;
+      SET_E_NTH_VERT(e,0,v_next);
       SET_E_NTH_TRI(e,0,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));
+      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);
 
-
-      if(qtinset(*delptr,t_next_id))
-	{
-#ifdef DEBUG
-	  eputs("smVertex_star_polygon(): id already in set\n");
-#endif	  
-	  free_list(elist);
-	  return(NULL);
-	}
-      else
-	qtaddelem(*delptr,t_next_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)
-SM *sm;
-FVECT v0,v1;
-LIST *l;
-{
-    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;
-
-    while(el)
-    {
-      e = (int)LIST_DATA(el);
-      id_e0 = E_NTH_VERT(e,0);
-      id_e1 = E_NTH_VERT(e,1);
-
-      VSUB(e0,SM_NTH_WV(sm,id_e0),SM_VIEW_CENTER(sm));
-      VSUB(e1,SM_NTH_WV(sm,id_e1),SM_VIEW_CENTER(sm));
-      if(sedge_intersect(v0,v1,e0,e1))
-	return(TRUE);
-
-      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 e,id;
-    LIST *el;
-    FVECT v;
-
-    /* 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 = 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(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 = INVALID;
-    list = lptr = *l;
-
-    if((new_e = eNew_edge())==INVALID)
-     {
-#ifdef DEBUG
-	    eputs("split_edge_list():Too many edges\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,INVALID);
-    SET_E_NTH_TRI(new_e,1,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);
-	}
-
-    }
-    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 = 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");
-#endif
-	  *l = NULL;
-	  return(FALSE);
-	}
-
-    }
-    end = lptr;
-    list = add_data_to_circular_list(list,&end,new_e);
-    *l = list;
-    return(TRUE);
-}
-
-
-int
-smTriangulate_convex(sm,plist,add_ptr)
-SM *sm;
-LIST *plist,**add_ptr;
-{
-    int t_id,e_id0,e_id1,e_id2;
-    int v_id0,v_id1,v_id2;
-    LIST *lptr;
-    int cnt;
-
-    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);
-	*add_ptr = push_data(*add_ptr,t_id);
-
-	/* 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,add_ptr)
-SM *sm;
-LIST *plist,**add_ptr;
-{
-    LIST *l,*el1;
-    FVECT v0,v1,v2;
-    int id0,id1,id2,e,id_next;
-    char flipped;
-    int done;
-
-    l = plist;
-    
-    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))
-      {
-	if(l == plist)
-	  break;
-	else
-	  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 == INVALID)
-      {
-#ifdef DEBUG
-	  eputs("smTriangulate_elist():Unable to find convex vertex\n");
-#endif
-	  return(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,add_ptr);
-      if(done)
-	done = smTriangulate_elist(sm,el1,add_ptr);
-      return(done);
-    }
-    done = smTriangulate_convex(sm,plist,add_ptr);
-    return(done);
-}
-
-int
 smTriangulate_add_tri(sm,id0,id1,id2,e0,e1,e2ptr)
 SM *sm;
 int id0,id1,id2,e0,e1,*e2ptr;
 {
-  int t_id;
-  int e2;
+  int t_id,e2;
+  TRI *t;
 
-  t_id = smAdd_tri(sm,id0,id1,id2);
+  t_id = smAdd_tri(sm,id0,id1,id2,&t);
   if(*e2ptr == 0)
   {
     e2 = eNew_edge();
@@ -448,145 +159,330 @@ int id0,id1,id2,e0,e1,*e2ptr;
   SET_E_NTH_TRI(e0,0,t_id);
   SET_E_NTH_TRI(e1,0,t_id);
   SET_E_NTH_TRI(e2,0,t_id);
-
+#ifdef DEBUG
+#if DEBUG > 1 
+  if(E_NTH_TRI(e0,1) == E_NTH_TRI(e0,0) ||
+     E_NTH_TRI(e1,1) == E_NTH_TRI(e1,0) ||
+     E_NTH_TRI(e2,1) == E_NTH_TRI(e2,0))
+    error(CONSISTENCY,"Non manifold\n");
+#endif
+#endif
   *e2ptr = e2;
   return(t_id);
+
 }
+
+int 
+smTriangulate_quad(sm,l)
+     SM *sm;
+     LIST *l;
+{
+  int e1,e2,e3,e4,e,t_id,id0,id1,id2,id3;
+  FVECT v0,v1,v2,v3,n;
+  double dp;
+  TRI *tri;
+  LIST *lptr;
+
+#ifdef DEBUG
+  if(LIST_NEXT(LIST_NEXT(LIST_NEXT(LIST_NEXT(l)))) != l)
+  {
+    eputs("smTriangulate_quad(): not a quadrilateral\n");
+    return(FALSE);
+  }
+    eputs("smTriangulate_quad():\n");
+#endif
+  lptr=l;
+  e1 = (int)LIST_DATA(l);
+  id0 = E_NTH_VERT(e1,0);
+  id1 = E_NTH_VERT(e1,1);
+  VSUB(v0,SM_NTH_WV(sm,id0),SM_VIEW_CENTER(sm));
+  VSUB(v1,SM_NTH_WV(sm,id1),SM_VIEW_CENTER(sm));   
+  /* Get v2 */
+  l = LIST_NEXT(l);
+  e2 = (int)LIST_DATA(l);
+  id2 = E_NTH_VERT(e2,1);
+  VSUB(v2,SM_NTH_WV(sm,id2),SM_VIEW_CENTER(sm));
+  /* Get v3 */
+  l = LIST_NEXT(l);
+  e3 = (int)LIST_DATA(l);
+  id3 = E_NTH_VERT(e3,1);
+  VSUB(v3,SM_NTH_WV(sm,id3),SM_VIEW_CENTER(sm));
+  l = LIST_NEXT(l);
+  e4 = (int)LIST_DATA(l);
+  free_list(lptr);
+
+  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);
+      }
+    }
+
+  }
+#ifdef DEBUG
+#if DEBUG > 1
+  VCROSS(n,v1,v3);
+  normalize(n);
+  normalize(v0);
+  normalize(v2);
+  dp = DOT(n,v2);
+  if((dp < 0) || (dp < EV_EPS) || (DOT(n,v0) >= 0.0))	
+    error(CONSISTENCY,"smTriangulate_quad: cannot triangulate\n");
+#endif
+#endif
+  t_id = smTriangulate_add_tri(sm,id1,id2,id3,e2,e3,&e);
+  e *= -1;
+  t_id = smTriangulate_add_tri(sm,id3,id0,id1,e4,e1,&e);
+  return(TRUE);
+}
+
+/* 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_elist_new(sm,id,plist,add_ptr)
+smTriangulate(sm,id,plist)
 SM *sm;
 int id;
-LIST *plist,**add_ptr;
+LIST *plist;
 {
     LIST *l,*prev,*t;
     FVECT v0,v1,v2,n,p;
-    int is_tri,loop,t_id,id0,id1,id2,e2,eprev,enext;
-    double dp;
+    int is_tri,cut,t_id,id0,id1,id2,e2,e1,enew;
+    double dp1,dp2;
 
-    smDir(sm,p,id);
-    enext=0;
-    is_tri= loop = FALSE;
-    l = prev = plist;
-    /* get v0,v1,v2 */
-    eprev = (int)LIST_DATA(l);
-    id0 = E_NTH_VERT(eprev,0);
-    id1 = E_NTH_VERT(eprev,1);
-    smDir(sm,v0,id0);
-    smDir(sm,v1,id1);   
+    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)
     {
       l = LIST_NEXT(l);
+      /* Get v2 */
       e2 = (int)LIST_DATA(l);
       id2 = E_NTH_VERT(e2,1);
-      /* Check if have a triangle */
+      VSUB(v2,SM_NTH_WV(sm,id2),SM_VIEW_CENTER(sm));
+      normalize(v2);
       if(LIST_NEXT(LIST_NEXT(l)) == prev)
-      {
-	is_tri = TRUE;
-	break;
+      {/* Check if have a triangle */
+	   is_tri = TRUE;	
+	   break;
       }
-      if(LIST_NEXT(l) == plist)
-      {
-	if(!loop)
-	  loop = 1;
-	else
-	  loop++;
-	if(loop > 3)
-	  break;
-      }
-      smDir(sm,v2,id2);
-      /* determine if convex (left turn), or concave(right turn) angle */
-      if(!convex_angle(v0,v1,v2))
-      {
-	VCOPY(v0,v1);
-	VCOPY(v1,v2);
-	id0 = id1;
-	id1 = id2;
-	prev = l;
-	eprev = e2;
-	continue;
-      }
+      /* determine if v0-v1-v2 is convex:defined clockwise on the sphere-
+       so switch orientation
+       */
       VCROSS(n,v0,v2);
-      dp = DOT(n,p);
-      if(loop <=1 && (!ZERO(dp) && dp  < 0.0))
+      normalize(n);
+      dp1 = DOT(n,p);      
+      if(dp1  <= 0.0)
       {
-	VCOPY(v0,v1);
-	VCOPY(v1,v2);
-	id0 = id1;
-	id1 = id2;
-	eprev = e2;
-	prev = l;
-	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;
+	}
       }
-      loop = FALSE;
-
-      enext = 0;
-      t_id = smTriangulate_add_tri(sm,id0,id1,id2,eprev,e2,&enext);
-      *add_ptr = push_data(*add_ptr,t_id);
-
-      LIST_NEXT(prev) = LIST_NEXT(l);
-      LIST_DATA(prev) = eprev = -enext;
-      LIST_NEXT(l)=NULL;
-      if(l== plist)
-	plist = prev;
-      free_list(l);
-      l = prev;
+      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)
+      {
+	if(LIST_NEXT(LIST_NEXT(l)) == prev)
+	  {/* Check if have a triangle */
+	    is_tri = TRUE;	
+	    break;
+	  }
+	if(!cut)
+	  break;
+        cut = FALSE;
+      }
+      prev = l;
     }
     if(is_tri)
     {
       l = LIST_NEXT(l);
-      enext = (int)LIST_DATA(l);
-      t_id = smTriangulate_add_tri(sm,id0,id1,id2,eprev,e2,&enext);
-      *add_ptr = push_data(*add_ptr,t_id);
+      enew = (int)LIST_DATA(l);
+      t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew);
       free_list(l);
-     }
+    }
     else
-      {
-#ifdef DEBUG       
-	eputs("smTriangulate_elist()Unable to triangulate\n");
+      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
-	return(FALSE);
-      }
-    return(TRUE);
-}
 
-int
-smTriangulate(sm,p_id,plist,add_ptr)
-SM *sm;
-int p_id;
-LIST *plist,**add_ptr;
-{
-    int e,id_t0,id_t1,e0,e1;
-    int test;
-
-    test = smTriangulate_elist_new(sm,p_id,plist,add_ptr);
+    /* Set triangle adjacencies based on edge adjacencies */ 
+    FOR_ALL_EDGES(enew)
+    {
+      id0 = E_NTH_TRI(enew,0);
+      id1 = E_NTH_TRI(enew,1);
+	
+      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
-    test = smTriangulate_elist(sm,plist,add_ptr);
+ {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
 
-    if(!test)
-       return(test);
+#ifdef DEBUG
+#if DEBUG > 1
+    {
+      TRI *t0,*t1,*n;
 
-    FOR_ALL_EDGES(e)
+    FOR_ALL_EDGES(enew)
     {
-	id_t0 = E_NTH_TRI(e,0);
-	id_t1 = E_NTH_TRI(e,1);
-	if((id_t0==INVALID) || (id_t1==INVALID))
-	{
-#ifdef DEBUG
-	   eputs("smTriangulate(): Unassigned edge neighbor\n");
-#endif
-	    continue;
-	}
-	
-	e0 = T_WHICH_V(SM_NTH_TRI(sm,id_t0),E_NTH_VERT(e,0));
-	T_NTH_NBR(SM_NTH_TRI(sm,id_t0),e0) = id_t1;
+      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");
 
-	e1 = T_WHICH_V(SM_NTH_TRI(sm,id_t1),E_NTH_VERT(e,1));
-	T_NTH_NBR(SM_NTH_TRI(sm,id_t1),e1) = id_t0;
+      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)
@@ -601,50 +497,193 @@ 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,add_list,delptr)
+
+void
+smTris_swap_edge(sm,t_id,t1_id,e,e1,tn_id,tn1_id)
    SM *sm;
-   LIST *add_list;
-   QUADTREE *delptr;
+   int t_id,t1_id;
+   int e,e1;
+   int *tn_id,*tn1_id;
+{
+    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;
+    
+    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)
     {
 	t0_id = E_NTH_TRI(e,0);
 	t1_id = E_NTH_TRI(e,1);
-	if((t0_id==INVALID) || (t1_id==INVALID))
-	{
 #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;
-	}
-	e0 = T_WHICH_V(SM_NTH_TRI(sm,t0_id),E_NTH_VERT(e,0));
-	e1 = T_WHICH_V(SM_NTH_TRI(sm,t1_id),E_NTH_VERT(-e,0));
-	e0_next = (e0+2)%3;
-	e1_next = (e1+2)%3;
+	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);
+
 	v0_id = E_NTH_VERT(e,0);
 	v1_id = E_NTH_VERT(e,1);
-	v2_id = T_NTH_V(SM_NTH_TRI(sm,t0_id),e0_next);
-	p_id = T_NTH_V(SM_NTH_TRI(sm,t1_id),e1_next);
+	v2_id = T_NTH_V(t0,e0);
+	p_id = T_NTH_V(t1,e1);
 
-	smDir_in_cone(sm,v0,v0_id);
-	smDir_in_cone(sm,v1,v1_id);
-	smDir_in_cone(sm,v2,v2_id);
+	smDir(sm,v0,v0_id);
+	smDir(sm,v1,v1_id);
+	smDir(sm,v2,v2_id);
 	
-	VCOPY(p,SM_NTH_WV(sm,p_id));	
-	VSUB(p,p,SM_VIEW_CENTER(sm));
-	if(point_in_cone(p,v0,v1,v2))
+	VSUB(p,SM_NTH_WV(sm,p_id),SM_VIEW_CENTER(sm));
+	normalize(p);
+	if(pt_in_cone(p,v2,v1,v0))
 	{
-	   smTris_swap_edge(sm,t0_id,t1_id,e0,e1,&t0_nid,&t1_nid,&add_list,
-			    delptr);
-	    
+	   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)==t0_id || E_NTH_TRI(i,0)==t1_id)
@@ -662,6 +701,10 @@ smFix_edges(sm,add_list,delptr)
 		else
 		  SET_E_NTH_TRI(i,1,t1_nid);
 	      }
+#ifdef DEBUG 
+	      if(E_NTH_TRI(i,1) == E_NTH_TRI(i,0) )
+		error(CONSISTENCY,"invalid edge\n");
+#endif
 	    }
 	    t0_id = t0_nid;
 	    t1_id = t1_nid;
@@ -670,57 +713,70 @@ smFix_edges(sm,add_list,delptr)
 	    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_list,qtqueryset(*delptr));
 }
 
+
+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 t_id;
-    LIST *elist,*add_list;
-    int cnt,debug;
-    QUADTREE delnode;
-    /* 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,&delnode);
-
-    if(!elist)
-    {
-#ifdef DEBUG
-      eputs("smMesh_remove_vertex(): Unable to remove vertex");
+    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
-      qtfreeleaf(delnode);
-      return(FALSE);
-    }
-    add_list = NULL;
-    /* Triangulate spherical polygon */
-    if(!smTriangulate(sm,id,elist,&add_list))
-    {
-      while(add_list)
-      {
-	t_id = pop_list(&add_list);
-	smDelete_tri(sm,t_id);
-      }
-      qtfreeleaf(delnode);
-      return(FALSE);
-    }
-    /* Fix up new triangles to be Delaunay */
-    smFix_edges(sm,add_list,&delnode);
 
-    qtfreeleaf(delnode);
+    /* Triangulate polygonal hole  */
+    smTriangulate(sm,id,b_list);
+    
+    /* Fix up new triangles to be Delaunay*/
+
+    smFixEdges(sm);
+    smDelete_samp(sm,id);
+    eClear_edges();
     return(TRUE);
 }
    
 
 
-
+