--- ray/src/hd/sm_del.c	1999/03/05 16:32:22	3.11
+++ ray/src/hd/sm_del.c	1999/06/10 15:22:22	3.12
@@ -11,34 +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;
 
-smFree_tri(sm,id)
+smFree_tri(sm,id,t)
 SM *sm;
 int id;
+TRI *t;
 {
-  TRI *tri,*t;
-
-  tri = SM_NTH_TRI(sm,id);
   /* Add to the free_list */
-  smClear_tri_flags(sm,id);
-  T_NEXT_FREE(tri) = SM_FREE_TRIS(sm);
+  T_NEXT_FREE(t) = SM_FREE_TRIS(sm);
   SM_FREE_TRIS(sm) = id;
-  T_VALID_FLAG(tri) = -1;
+#ifdef DEBUG
+  T_VALID_FLAG(t) = INVALID;
+#endif
 }
 
 /* Assumes mesh pointers have been cleaned up appropriately: just deletes from
    Point location and triangle data structure
 */
-smDelete_tri(sm,t_id)
+smDelete_tri(sm,t_id,t)
 SM *sm;
 int t_id;
+TRI *t;
 {
 
   /* NOTE: Assumes that a new triangle adjacent to each vertex
@@ -47,16 +46,11 @@ int t_id;
      to id because the vertices can no longer
      point to tri id as being the first triangle pointer
   */
-  SM_SAMPLE_TRIS(sm)--;
-  if(SM_IS_NTH_T_NEW(sm,t_id))
-     smNew_tri_cnt--;
+  SM_CLR_NTH_T_ACTIVE(sm,t_id);
+  smFree_tri(sm,t_id,t);
 
-  smClear_tri_flags(sm,t_id);
-
-  smFree_tri(sm,t_id);
 }
 
-
 int
 eNew_edge()
 {
@@ -66,35 +60,36 @@ eNew_edge()
 
   if(Ecnt >= Max_edges)
     {
+#ifdef DEBUG
       if(Max_edges > 10000)
-      {
-	eputs("Too many edges in vertex loop\n");
-	return(-1);
-      }
+	error(CONSISTENCY,"Too many edges in vertex loop\n");
+#endif
       Max_edges += 100;
       if(!(Edges = (EDGE *)realloc(Edges,(Max_edges+1)*sizeof(EDGE))))
 	goto memerr;
     }
+#ifdef DEBUG
+  SET_E_NTH_TRI(Ecnt+1,0,INVALID);
+  SET_E_NTH_TRI(Ecnt+1,1,INVALID);
+#endif
   return(++Ecnt);
 
 memerr:
-  error(SYSTEM,"eNew_edge(): Unable to allocate memory");
+  error(SYSTEM,"eNew_edge): Unable to allocate memory");
 }
 
 /* Return list of edges defining polygon formed by boundary of triangles 
 adjacent to id. Return set of triangles adjacent to id to delete in delptr
 */
 LIST 
-*smVertexStar(sm,id,del_ptr)
+*smVertexStar(sm,id)
 SM *sm;
 int id;
-LIST **del_ptr;
 {
     TRI *tri,*t_next;
     LIST *elist,*end;
-    int e,t_id,v_next,t_next_id,b_id,v_id;
+    int e,t_id,v_next,t_next_id,b_id,v_id,t_last_id;
 
-    eClear_edges();
     elist = end =  NULL;
 
     /* Get the first triangle adjacent to vertex id */
@@ -102,7 +97,6 @@ LIST **del_ptr;
     tri = SM_NTH_TRI(sm,t_id);
 
     e = eNew_edge();
-
     /* Get the  next vertex on the polygon boundary */
     v_id = T_WHICH_V(tri,id);
     b_id = (v_id + 1)%3;
@@ -110,21 +104,25 @@ LIST **del_ptr;
     SET_E_NTH_VERT(e,0,T_NTH_V(tri,b_id));
     SET_E_NTH_TRI(e,0,INVALID);
     SET_E_NTH_TRI(e,1,T_NTH_NBR(tri,v_id));
-    v_next = T_NTH_V(tri,(b_id+1)%3);
+     v_next = T_NTH_V(tri,(b_id+1)%3);
     SET_E_NTH_VERT(e,1,v_next);
+  
     elist = add_data_to_circular_list(elist,&end,e);
     t_next_id = t_id;
     t_next = tri;
-
-    *del_ptr = push_data(*del_ptr,t_id);
+    t_last_id = -1;
+ 
     /* Create a set to hold all of the triangles for deletion later */
 
     while((t_next_id = T_NTH_NBR(t_next,b_id)) != t_id)
     {
       e = eNew_edge();
-      if(e== INVALID)
-	return(NULL);
+      if(t_last_id != -1)
+      {
+	smDelete_tri(sm,t_last_id,t_next);
+      }
       t_next = SM_NTH_TRI(sm,t_next_id);
+      t_last_id = t_next_id;
       SET_E_NTH_VERT(e,0,v_next);
       SET_E_NTH_TRI(e,0,INVALID);
       v_id = T_WHICH_V(t_next,id);
@@ -133,25 +131,22 @@ LIST **del_ptr;
       v_next = T_NTH_V(t_next,(b_id+1)%3);
       SET_E_NTH_VERT(e,1,v_next);
       elist = add_data_to_circular_list(elist,&end,e);
-      *del_ptr = push_data(*del_ptr,t_next_id);
+
     }
-    return(elist);
+    smDelete_tri(sm,t_last_id,t_next);
+    smDelete_tri(sm,t_id,tri); 
+   return(elist);
 }
 
-
 int
 smTriangulate_add_tri(sm,id0,id1,id2,e0,e1,e2ptr)
 SM *sm;
 int id0,id1,id2,e0,e1,*e2ptr;
 {
-  int t_id;
-  int e2;
+  int t_id,e2;
+  TRI *t;
 
-#ifdef DEBUG 
-  if(id0 == INVALID || id1==INVALID || id2==INVALID)
-    error(CONSISTENCY,"bad id- smTriangulate_add_tri()\n");
-#endif
-  t_id = smAdd_tri(sm,id0,id1,id2);
+  t_id = smAdd_tri(sm,id0,id1,id2,&t);
   if(*e2ptr == 0)
   {
     e2 = eNew_edge();
@@ -164,48 +159,95 @@ int id0,id1,id2,e0,e1,*e2ptr;
   SET_E_NTH_TRI(e0,0,t_id);
   SET_E_NTH_TRI(e1,0,t_id);
   SET_E_NTH_TRI(e2,0,t_id);
-
+#ifdef DEBUG
+#if DEBUG > 1 
+  if(E_NTH_TRI(e0,1) == E_NTH_TRI(e0,0) ||
+     E_NTH_TRI(e1,1) == E_NTH_TRI(e1,0) ||
+     E_NTH_TRI(e2,1) == E_NTH_TRI(e2,0))
+    error(CONSISTENCY,"Non manifold\n");
+#endif
+#endif
   *e2ptr = e2;
   return(t_id);
+
 }
 
-int
-smTriangulateConvex(sm,plist,add_ptr)
-SM *sm;
-LIST *plist,**add_ptr;
+int 
+smTriangulate_quad(sm,l)
+     SM *sm;
+     LIST *l;
 {
-    int t_id,e_id0,e_id1,e_id2;
-    int v_id0,v_id1,v_id2;
-    LIST *lptr;
+  int e1,e2,e3,e4,e,t_id,id0,id1,id2,id3;
+  FVECT v0,v1,v2,v3,n;
+  double dp;
+  TRI *tri;
+  LIST *lptr;
 
-    lptr = plist;
-    e_id0 = (int)LIST_DATA(lptr);
-    v_id0 = E_NTH_VERT(e_id0,0);
-    lptr = LIST_NEXT(lptr);
-    while(LIST_NEXT(lptr) != plist)
-    {
-	e_id1 = (int)LIST_DATA(lptr);
-	v_id1 = E_NTH_VERT(e_id1,0);
-	v_id2 = E_NTH_VERT(e_id1,1);
-	lptr = LIST_NEXT(lptr);
+#ifdef DEBUG
+  if(LIST_NEXT(LIST_NEXT(LIST_NEXT(LIST_NEXT(l)))) != l)
+  {
+    eputs("smTriangulate_quad(): not a quadrilateral\n");
+    return(FALSE);
+  }
+    eputs("smTriangulate_quad():\n");
+#endif
+  lptr=l;
+  e1 = (int)LIST_DATA(l);
+  id0 = E_NTH_VERT(e1,0);
+  id1 = E_NTH_VERT(e1,1);
+  VSUB(v0,SM_NTH_WV(sm,id0),SM_VIEW_CENTER(sm));
+  VSUB(v1,SM_NTH_WV(sm,id1),SM_VIEW_CENTER(sm));   
+  /* Get v2 */
+  l = LIST_NEXT(l);
+  e2 = (int)LIST_DATA(l);
+  id2 = E_NTH_VERT(e2,1);
+  VSUB(v2,SM_NTH_WV(sm,id2),SM_VIEW_CENTER(sm));
+  /* Get v3 */
+  l = LIST_NEXT(l);
+  e3 = (int)LIST_DATA(l);
+  id3 = E_NTH_VERT(e3,1);
+  VSUB(v3,SM_NTH_WV(sm,id3),SM_VIEW_CENTER(sm));
+  l = LIST_NEXT(l);
+  e4 = (int)LIST_DATA(l);
+  free_list(lptr);
 
-	if(LIST_NEXT(lptr) != plist)	
-	  e_id2 = 0;
-	else
-	   e_id2 = (int)LIST_DATA(lptr);
-	t_id = smTriangulate_add_tri(sm,v_id0,v_id1,v_id2,e_id0,e_id1,&e_id2);
-	*add_ptr = push_data(*add_ptr,t_id);
-	e_id0 = -e_id2;
+  VCROSS(n,v0,v2);
+  normalize(n);
+  normalize(v1);
+  dp = DOT(n,v1);
+  e = 0;
+  if(dp > 0) 
+  {
+    if(dp  >= EV_EPS)
+    {
+      normalize(v3);
+      if(DOT(n,v3) < 0)	
+      {
+	t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&e);
+	e *= -1;
+	t_id = smTriangulate_add_tri(sm,id2,id3,id0,e3,e4,&e);
+	return(TRUE);
+      }
     }
-    free_list(plist);
-    return(TRUE);
-}
-#ifdef TEST_DRIVER
-FVECT Norm[500],B_V[500];
-int Ncnt,Bcnt,Del=0;
+
+  }
+#ifdef DEBUG
+#if DEBUG > 1
+  VCROSS(n,v1,v3);
+  normalize(n);
+  normalize(v0);
+  normalize(v2);
+  dp = DOT(n,v2);
+  if((dp < 0) || (dp < EV_EPS) || (DOT(n,v0) >= 0.0))	
+    error(CONSISTENCY,"smTriangulate_quad: cannot triangulate\n");
 #endif
+#endif
+  t_id = smTriangulate_add_tri(sm,id1,id2,id3,e2,e3,&e);
+  e *= -1;
+  t_id = smTriangulate_add_tri(sm,id3,id0,id1,e4,e1,&e);
+  return(TRUE);
+}
 
-
 /* Triangulate the polygon defined by plist, and generating vertex p_id.
    Return list of added triangles in list add_ptr. Returns TRUE if 
    successful, FALSE otherwise. This is NOT a general triangulation routine,
@@ -213,36 +255,31 @@ int Ncnt,Bcnt,Del=0;
 */
 
 int
-smTriangulate(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,is_convex,cut,t_id,id0,id1,id2,e2,e1,enew;
-    double dp;
-    static int debug=0;
+    int is_tri,cut,t_id,id0,id1,id2,e2,e1,enew;
+    double dp1,dp2;
 
     VSUB(p,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
-    enew = 0;
-    is_convex = TRUE;
-    cut = is_tri= FALSE;
-    l = prev = plist;
+    normalize(p);
+    l = plist;
 
+    enew = 0;
+    cut = is_tri=  FALSE;
+    prev = l;
     /* get v0,v1 */
     e1 = (int)LIST_DATA(l);
     id0 = E_NTH_VERT(e1,0);
     id1 = E_NTH_VERT(e1,1);
     VSUB(v0,SM_NTH_WV(sm,id0),SM_VIEW_CENTER(sm));
+    normalize(v0);
     VSUB(v1,SM_NTH_WV(sm,id1),SM_VIEW_CENTER(sm));   
-#ifdef TEST_DRIVER
-    Del = TRUE;
-    VCOPY(B_V[0],v0);
-    VCOPY(B_V[1],v1);
-    Bcnt = 2;
-    Ncnt = 0;
-#endif
+    normalize(v1);
     while(l)
     {
       l = LIST_NEXT(l);
@@ -250,50 +287,48 @@ LIST *plist,**add_ptr;
       e2 = (int)LIST_DATA(l);
       id2 = E_NTH_VERT(e2,1);
       VSUB(v2,SM_NTH_WV(sm,id2),SM_VIEW_CENTER(sm));
-#ifdef TEST_DRIVER
-      VCOPY(B_V[Bcnt++],v2);
-#endif
+      normalize(v2);
       if(LIST_NEXT(LIST_NEXT(l)) == prev)
       {/* Check if have a triangle */
 	   is_tri = TRUE;	
 	   break;
       }
-
       /* determine if v0-v1-v2 is convex:defined clockwise on the sphere-
        so switch orientation
        */
-      if(convex_angle(v2,v1,v0))
+      VCROSS(n,v0,v2);
+      normalize(n);
+      dp1 = DOT(n,p);      
+      if(dp1  <= 0.0)
       {
 	  /* test if safe to cut off v0-v1-v2 by testing if p lies outside of
          triangle v0-v1-v2: if so, because plist is the star polygon around p,
           the new edge v2-v0 cannot intersect any existing edges
         */
-	VCROSS(n,v0,v2);
-        dp = DOT(n,p);
-        if(dp  <= 0.0)
-	{
+	dp1 = DOT(n,v1);
+	/* Distance from point s to plane is d = |N.p0s|/||N|| */
+        /* sp1 = s-p0 for point on plane p0, a.(b+c) = a.b + a.c */
+	if(dp1 >= EV_EPS)
+        {
 	    /* remove edges e1,e2 and add triangle id0,id1,id2 */
-          enew = 0;
-          t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew);
-	  cut = TRUE;
-          *add_ptr = push_data(*add_ptr,t_id);
-           /* Insert edge enew into the list, reuse prev list element */
-	  LIST_NEXT(prev) = LIST_NEXT(l);
-	  LIST_DATA(prev) = e1 = -enew;
-           /* If removing head of list- reset plist pointer */
-          if(l== plist)
-	    plist = prev;
-	  /* free list element for e2 */
-	  LIST_NEXT(l)=NULL;
-	  free_list(l);
-	  l = prev;
-	  VCOPY(v1,v2);
-	  id1 = id2;
-	  continue;
+	    enew = 0;
+	    t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew);
+	    cut = TRUE;
+	    /* Insert edge enew into the list, reuse prev list element */
+	    LIST_NEXT(prev) = LIST_NEXT(l);
+	    LIST_DATA(prev) = e1 = -enew;
+	    /* If removing head of list- reset plist pointer */
+	    if(l== plist)
+	      plist = prev;
+	    /* free list element for e2 */
+	    LIST_NEXT(l)=NULL;
+	    free_list(l);
+	    l = prev;
+	    VCOPY(v1,v2);
+	    id1 = id2;
+	    continue;
 	}
       }
-      else
-	is_convex = FALSE;
       VCOPY(v0,v1);
       VCOPY(v1,v2);
       id0 = id1;
@@ -308,25 +343,9 @@ LIST *plist,**add_ptr;
 	    is_tri = TRUE;	
 	    break;
 	  }
-
-	if(is_convex)
-	  break;
 	if(!cut)
-	{
-#ifdef DEBUG
-	  eputs("smTriangulate():Unable to triangulate\n");
-#endif
-	 free_list(l);
-	  while(*add_ptr)
-	  {
-	    t_id = pop_list(add_ptr);
-	    smDelete_tri(sm,t_id);
-	  }
-	  return(FALSE);
-	 }
-        
+	  break;
         cut = FALSE;
-	is_convex = TRUE;
       }
       prev = l;
     }
@@ -335,12 +354,20 @@ LIST *plist,**add_ptr;
       l = LIST_NEXT(l);
       enew = (int)LIST_DATA(l);
       t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew);
-      *add_ptr = push_data(*add_ptr,t_id);
       free_list(l);
     }
     else
-      if(!smTriangulateConvex(sm,l,add_ptr))
-	return(FALSE);
+      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
 
     /* Set triangle adjacencies based on edge adjacencies */ 
     FOR_ALL_EDGES(enew)
@@ -354,7 +381,108 @@ LIST *plist,**add_ptr;
       e2 = (T_WHICH_V(SM_NTH_TRI(sm,id1),E_NTH_VERT(enew,1))+2)%3;
       T_NTH_NBR(SM_NTH_TRI(sm,id1),e2) = id0;
     }
+#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
+
+#ifdef DEBUG
+#if DEBUG > 1
+    {
+      TRI *t0,*t1,*n;
+
+    FOR_ALL_EDGES(enew)
+    {
+      id0 = E_NTH_TRI(enew,0);
+      id1 = E_NTH_TRI(enew,1);
+      t0 = SM_NTH_TRI(sm,id0);
+      t1 = SM_NTH_TRI(sm,id1);
+      if(T_NTH_NBR(t0,0) == T_NTH_NBR(t0,1)  || 
+       T_NTH_NBR(t0,1) == T_NTH_NBR(t0,2)  ||
+       T_NTH_NBR(t0,0) == T_NTH_NBR(t0,2)) 
+      error(CONSISTENCY,"Invalid topology\n");
+    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,0))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,1))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,2))))
+      error(CONSISTENCY,"Invalid topology\n");
+    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,0));
+    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  || 
+       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
+       T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) 
+      error(CONSISTENCY,"Invalid topology\n");
+    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
+      error(CONSISTENCY,"Invalid topology\n");
+    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,1));
+    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  || 
+       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
+       T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) 
+      error(CONSISTENCY,"Invalid topology\n");
+    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
+      error(CONSISTENCY,"Invalid topology\n");
+    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,2));
+    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  || 
+       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
+       T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) 
+      error(CONSISTENCY,"Invalid topology\n");
+    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
+      error(CONSISTENCY,"Invalid topology\n");
+
+      if(T_NTH_NBR(t1,0) == T_NTH_NBR(t1,1)  || 
+       T_NTH_NBR(t1,1) == T_NTH_NBR(t1,2)  ||
+       T_NTH_NBR(t1,0) == T_NTH_NBR(t1,2)) 
+      error(CONSISTENCY,"Invalid topology\n");
+    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,0))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,1))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,2))))
+      error(CONSISTENCY,"Invalid topology\n");
+    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,0));
+    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  || 
+       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
+       T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) 
+      error(CONSISTENCY,"Invalid topology\n");
+    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
+      error(CONSISTENCY,"Invalid topology\n");
+    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,1));
+    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  || 
+       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
+       T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) 
+      error(CONSISTENCY,"Invalid topology\n");
+    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
+      error(CONSISTENCY,"Invalid topology\n");
+    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,2));
+    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  || 
+       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
+       T_NTH_NBR(n,0) == T_NTH_NBR(n,2)) 
+      error(CONSISTENCY,"Invalid topology\n");
+    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
+       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
+      error(CONSISTENCY,"Invalid topology\n");
+    }
+    }
+#endif
+#endif
     return(TRUE);
+
+Terror:	  
+#ifdef DEBUG
+	  error(CONSISTENCY,"smTriangulate():Unable to triangulate\n");
+#endif
 }
 
 eIn_tri(e,t)
@@ -373,15 +501,153 @@ TRI *t;
   return(FALSE);
 }
 
+
+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 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,add_list)
+smFixEdges(sm)
    SM *sm;
-   LIST *add_list;
 {
     int e,t0_id,t1_id,e_new,e0,e1,e0_next,e1_next;
     int i,v0_id,v1_id,v2_id,p_id,t0_nid,t1_nid;
@@ -392,12 +658,10 @@ smFixEdges(sm,add_list)
     {
 	t0_id = E_NTH_TRI(e,0);
 	t1_id = E_NTH_TRI(e,1);
-	if((t0_id==INVALID) || (t1_id==INVALID))
-	{
 #ifdef DEBUG
-	    error(CONSISTENCY,"smFix_edges: Unassigned edge nbr\n");
+	if((t0_id==INVALID) || (t1_id==INVALID))
+	  error(CONSISTENCY,"smFix_edges: Unassigned edge nbr\n");
 #endif
-	}
 	t0 = SM_NTH_TRI(sm,t0_id);
 	t1 = SM_NTH_TRI(sm,t1_id);
 	e0 = T_NTH_NBR_PTR(t1_id,t0);
@@ -412,12 +676,11 @@ smFixEdges(sm,add_list)
 	smDir(sm,v1,v1_id);
 	smDir(sm,v2,v2_id);
 	
-	VCOPY(p,SM_NTH_WV(sm,p_id));	
-	VSUB(p,p,SM_VIEW_CENTER(sm));
-	if(point_in_cone(p,v0,v1,v2))
+	VSUB(p,SM_NTH_WV(sm,p_id),SM_VIEW_CENTER(sm));
+	normalize(p);
+	if(pt_in_cone(p,v2,v1,v0))
 	{
-	   smTris_swap_edge(sm,t0_id,t1_id,e0,e1,&t0_nid,&t1_nid,&add_list);
-	    
+	   smTris_swap_edge(sm,t0_id,t1_id,e0,e1,&t0_nid,&t1_nid);
 	   /* Adjust the triangle pointers of the remaining edges to be
 	      processed
             */
@@ -438,6 +701,10 @@ smFixEdges(sm,add_list)
 		else
 		  SET_E_NTH_TRI(i,1,t1_nid);
 	      }
+#ifdef DEBUG 
+	      if(E_NTH_TRI(i,1) == E_NTH_TRI(i,0) )
+		error(CONSISTENCY,"invalid edge\n");
+#endif
 	    }
 	    t0_id = t0_nid;
 	    t1_id = t1_nid;
@@ -446,12 +713,34 @@ smFixEdges(sm,add_list)
 	    SET_E_NTH_VERT(e_new,1,v2_id);
 	    SET_E_NTH_TRI(e_new,0,t0_id);
 	    SET_E_NTH_TRI(e_new,1,t1_id);
+#ifdef DEBUG 
+	      if(E_NTH_TRI(i,1) == E_NTH_TRI(i,0) )
+		error(CONSISTENCY,"invalid edge\n");
+#endif
 	}
     }
-    /* Add/Delete the appropriate triangles from the stree */
-    smUpdate_locator(sm,add_list);
 }
 
+
+smDelete_samp(sm,s_id)
+SM *sm;
+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.
  */
@@ -460,41 +749,28 @@ smRemoveVertex(sm,id)
    SM *sm;
    int id;
 {
-    LIST *b_list,*add_list,*del_list;
-    int t_id,i;
-    static int cnt=0;
-    OBJECT *optr,*os;
+    LIST *b_list;
     /* generate list of edges that form the boundary of the
-       polygon formed by the triangles adjacent to vertex 'id'
-     */
-    del_list = NULL;
-    b_list = smVertexStar(sm,id,&del_list);
-    if(!b_list)
-    {
-      if(del_list)
-	free_list(del_list);
-      return(FALSE);
-    }
-    add_list = NULL;
+       polygon formed by the triangles adjacent to vertex 'id'*/
+    b_list = smVertexStar(sm,id);
+#if 0
+ {int i;
+    eputs("\n\n");
+    for(i=1;i<=Ecnt;i++)
+      fprintf(stderr,"%d verts %d %d tris  %d %d\n",
+	      i,Edges[i].verts[0],Edges[i].verts[1],
+	      Edges[i].tris[0],Edges[i].tris[1]);
+ }
+#endif
+
     /* Triangulate polygonal hole  */
-    if(!smTriangulate(sm,id,b_list,&add_list))
-    {
-      free_list(del_list);
-      return(FALSE);
-    }
-    else
-    {
-      while(del_list)
-      {
-	t_id = pop_list(&del_list);
-	smDelete_tri(sm,t_id);
-      }     
-    }      
-    /* Fix up new triangles to be Delaunay-delnode contains set of
-     triangles to delete,add_list is the set of new triangles to add
-     */
-    smFixEdges(sm,add_list);
+    smTriangulate(sm,id,b_list);
+    
+    /* Fix up new triangles to be Delaunay*/
 
+    smFixEdges(sm);
+    smDelete_samp(sm,id);
+    eClear_edges();
     return(TRUE);
 }