src/hd/sm_del.c

/* Copyright (c) 1998 Silicon Graphics, Inc. */

#ifndef lint
static char SCCSid[] = "$SunId$ SGI";
#endif

/*
 *  sm_del.c
 */
#include "standard.h"
#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;
static EDGE *Edges=NULL;
static int Ecnt=0;

smFree_tri(sm,id)
SM *sm;
int id;
{
  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);
  SM_FREE_TRIS(sm) = id;
  T_VALID_FLAG(tri) = -1;
}

/* Assumes mesh pointers have been cleaned up appropriately: just deletes from
   Point location and triangle data structure
*/
smDelete_tri(sm,t_id)
SM *sm;
int t_id;
{

  /* 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
  */
  SM_SAMPLE_TRIS(sm)--;
  if(SM_IS_NTH_T_NEW(sm,t_id))
     smNew_tri_cnt--;

  smClear_tri_flags(sm,t_id);

  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)
    {
      if(Max_edges > 10000)
      {
        eputs("Too many edges in vertex loop\n");
        return(-1);
      }
      Max_edges += 100;
      if(!(Edges = (EDGE *)realloc(Edges,(Max_edges+1)*sizeof(EDGE))))
        goto memerr;
    }
  return(++Ecnt);

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,del_ptr)
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;

    eClear_edges();
    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);

    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);
    t_next_id = t_id;
    t_next = tri;

    *del_ptr = push_data(*del_ptr,t_id);
    /* 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);
      t_next = SM_NTH_TRI(sm,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);
      *del_ptr = push_data(*del_ptr,t_next_id);
    }
    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;

#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);
  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);

  *e2ptr = e2;
  return(t_id);
}

int
smTriangulateConvex(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;

    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);

        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;
    }
    free_list(plist);
    return(TRUE);
}
#ifdef TEST_DRIVER
FVECT Norm[500],B_V[500];
int Ncnt,Bcnt,Del=0;
#endif


/* 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(sm,id,plist,add_ptr)
SM *sm;
int id;
LIST *plist,**add_ptr;
{
    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;

    VSUB(p,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
    enew = 0;
    is_convex = TRUE;
    cut = is_tri= FALSE;
    l = prev = plist;

    /* 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));
    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
    while(l)
    {
      l = LIST_NEXT(l);
      /* Get v2 */
      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
      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))
      {
          /* 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)
        {
            /* 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;
        }
      }
      else
        is_convex = FALSE;
      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(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);
         }
        
        cut = FALSE;
        is_convex = TRUE;
      }
      prev = l;
    }
    if(is_tri)
    {
      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);

    /* 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;
    }
    return(TRUE);
}

eIn_tri(e,t)
int e;
TRI *t;
{

  if(T_NTH_V(t,0)==E_NTH_VERT(e,0)) 
    return(T_NTH_V(t,1)==E_NTH_VERT(e,1)||T_NTH_V(t,2)==E_NTH_VERT(e,1));
  else
    if(T_NTH_V(t,1)==E_NTH_VERT(e,0))
      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);
}

/* 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)
   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;
    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
            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);
        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(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,p_id));    
        VSUB(p,p,SM_VIEW_CENTER(sm));
        if(point_in_cone(p,v0,v1,v2))
        {
           smTris_swap_edge(sm,t0_id,t1_id,e0,e1,&t0_nid,&t1_nid,&add_list);
            
           /* 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)
              {
                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,t1_nid);
              }

              if(E_NTH_TRI(i,1)==t0_id || E_NTH_TRI(i,1)==t1_id)
              {
                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,t1_nid);
              }
            }
            t0_id = t0_nid;
            t1_id = t1_nid;
            e_new = eNew_edge();
            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);
        }
    }
    /* Add/Delete the appropriate triangles from the stree */
    smUpdate_locator(sm,add_list);
}

/* Remove vertex "id" from the mesh- and retriangulate the resulting
   hole: Returns TRUE if successful, FALSE otherwise.
 */
int
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;
    /* 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;
    /* 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);

    return(TRUE);
}
   

Revision:	3.11
Committed:	Fri Mar 5 16:32:22 1999 UTC (26 years, 7 months ago) by gwlarson
Content type:	text/plain
Branch:	MAIN
Changes since 3.10:	+5 -5 lines
Log Message:	small clean-up changes
#	User	Rev	Content
1	gwlarson	3.1	/* Copyright (c) 1998 Silicon Graphics, Inc. */
2
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ SGI";
5			#endif
6
7			/*
8			* sm_del.c
9			*/
10			#include "standard.h"
11	gwlarson	3.6	#include "sm_flag.h"
12	gwlarson	3.1	#include "sm_list.h"
13			#include "sm_geom.h"
14	gwlarson	3.6	#include "sm_qtree.h"
15			#include "sm_stree.h"
16	gwlarson	3.1	#include "sm.h"
17
18	gwlarson	3.6	static int Max_edges=200;
19			static EDGE *Edges=NULL;
20	gwlarson	3.1	static int Ecnt=0;
21
22			smFree_tri(sm,id)
23			SM *sm;
24			int id;
25			{
26	gwlarson	3.8	TRI tri,t;
27	gwlarson	3.1
28			tri = SM_NTH_TRI(sm,id);
29			/* Add to the free_list */
30	gwlarson	3.10	smClear_tri_flags(sm,id);
31	gwlarson	3.1	T_NEXT_FREE(tri) = SM_FREE_TRIS(sm);
32			SM_FREE_TRIS(sm) = id;
33			T_VALID_FLAG(tri) = -1;
34			}
35
36			/* Assumes mesh pointers have been cleaned up appropriately: just deletes from
37			Point location and triangle data structure
38			*/
39			smDelete_tri(sm,t_id)
40			SM *sm;
41			int t_id;
42			{
43
44			/* NOTE: Assumes that a new triangle adjacent to each vertex
45			has been added- before the deletion: replacing
46			the old tri- and therefore dont need to dereference any pointers
47			to id because the vertices can no longer
48			point to tri id as being the first triangle pointer
49			*/
50	gwlarson	3.9	SM_SAMPLE_TRIS(sm)--;
51	gwlarson	3.10	if(SM_IS_NTH_T_NEW(sm,t_id))
52			smNew_tri_cnt--;
53	gwlarson	3.9
54	gwlarson	3.1	smClear_tri_flags(sm,t_id);
55
56			smFree_tri(sm,t_id);
57	gwlarson	3.6	}
58	gwlarson	3.1
59	gwlarson	3.6
60			int
61			eNew_edge()
62			{
63			if(!Edges)
64			if(!(Edges = (EDGE )realloc(NULL,(Max_edges+1)sizeof(EDGE))))
65			goto memerr;
66
67			if(Ecnt >= Max_edges)
68			{
69			if(Max_edges > 10000)
70	gwlarson	3.10	{
71			eputs("Too many edges in vertex loop\n");
72			return(-1);
73			}
74	gwlarson	3.6	Max_edges += 100;
75			if(!(Edges = (EDGE )realloc(Edges,(Max_edges+1)sizeof(EDGE))))
76			goto memerr;
77			}
78			return(++Ecnt);
79
80			memerr:
81			error(SYSTEM,"eNew_edge(): Unable to allocate memory");
82	gwlarson	3.1	}
83
84	gwlarson	3.7	/* Return list of edges defining polygon formed by boundary of triangles
85			adjacent to id. Return set of triangles adjacent to id to delete in delptr
86			*/
87	gwlarson	3.5	LIST
88	gwlarson	3.11	*smVertexStar(sm,id,del_ptr)
89	gwlarson	3.1	SM *sm;
90			int id;
91	gwlarson	3.8	LIST **del_ptr;
92	gwlarson	3.1	{
93			TRI tri,t_next;
94	gwlarson	3.5	LIST elist,end;
95	gwlarson	3.7	int e,t_id,v_next,t_next_id,b_id,v_id;
96	gwlarson	3.1
97	gwlarson	3.7	eClear_edges();
98	gwlarson	3.1	elist = end = NULL;
99	gwlarson	3.7
100	gwlarson	3.1	/* Get the first triangle adjacent to vertex id */
101			t_id = SM_NTH_VERT(sm,id);
102			tri = SM_NTH_TRI(sm,t_id);
103
104	gwlarson	3.7	e = eNew_edge();
105	gwlarson	3.10
106	gwlarson	3.7	/* Get the next vertex on the polygon boundary */
107			v_id = T_WHICH_V(tri,id);
108			b_id = (v_id + 1)%3;
109			/* Create an edge */
110			SET_E_NTH_VERT(e,0,T_NTH_V(tri,b_id));
111	gwlarson	3.6	SET_E_NTH_TRI(e,0,INVALID);
112	gwlarson	3.7	SET_E_NTH_TRI(e,1,T_NTH_NBR(tri,v_id));
113			v_next = T_NTH_V(tri,(b_id+1)%3);
114			SET_E_NTH_VERT(e,1,v_next);
115	gwlarson	3.6	elist = add_data_to_circular_list(elist,&end,e);
116	gwlarson	3.1	t_next_id = t_id;
117			t_next = tri;
118
119	gwlarson	3.8	del_ptr = push_data(del_ptr,t_id);
120	gwlarson	3.7	/* Create a set to hold all of the triangles for deletion later */
121	gwlarson	3.1
122	gwlarson	3.7	while((t_next_id = T_NTH_NBR(t_next,b_id)) != t_id)
123			{
124			e = eNew_edge();
125	gwlarson	3.10	if(e== INVALID)
126			return(NULL);
127	gwlarson	3.6	t_next = SM_NTH_TRI(sm,t_next_id);
128	gwlarson	3.7	SET_E_NTH_VERT(e,0,v_next);
129	gwlarson	3.6	SET_E_NTH_TRI(e,0,INVALID);
130	gwlarson	3.7	v_id = T_WHICH_V(t_next,id);
131			b_id = (v_id + 1)%3;
132			SET_E_NTH_TRI(e,1,T_NTH_NBR(t_next,v_id));
133			v_next = T_NTH_V(t_next,(b_id+1)%3);
134			SET_E_NTH_VERT(e,1,v_next);
135	gwlarson	3.6	elist = add_data_to_circular_list(elist,&end,e);
136	gwlarson	3.8	del_ptr = push_data(del_ptr,t_next_id);
137	gwlarson	3.1	}
138			return(elist);
139			}
140
141	gwlarson	3.7
142	gwlarson	3.1	int
143	gwlarson	3.7	smTriangulate_add_tri(sm,id0,id1,id2,e0,e1,e2ptr)
144	gwlarson	3.1	SM *sm;
145	gwlarson	3.7	int id0,id1,id2,e0,e1,*e2ptr;
146	gwlarson	3.1	{
147	gwlarson	3.7	int t_id;
148			int e2;
149	gwlarson	3.1
150	gwlarson	3.7	#ifdef DEBUG
151			if(id0 == INVALID \|\| id1==INVALID \|\| id2==INVALID)
152			error(CONSISTENCY,"bad id- smTriangulate_add_tri()\n");
153	gwlarson	3.1	#endif
154	gwlarson	3.7	t_id = smAdd_tri(sm,id0,id1,id2);
155			if(*e2ptr == 0)
156			{
157			e2 = eNew_edge();
158			SET_E_NTH_VERT(e2,0,id2);
159			SET_E_NTH_VERT(e2,1,id0);
160			}
161			else
162			e2 = *e2ptr;
163			/* set appropriate tri for each edge*/
164			SET_E_NTH_TRI(e0,0,t_id);
165			SET_E_NTH_TRI(e1,0,t_id);
166			SET_E_NTH_TRI(e2,0,t_id);
167	gwlarson	3.1
168	gwlarson	3.7	*e2ptr = e2;
169			return(t_id);
170	gwlarson	3.1	}
171
172	gwlarson	3.3	int
173	gwlarson	3.7	smTriangulateConvex(sm,plist,add_ptr)
174	gwlarson	3.1	SM *sm;
175	gwlarson	3.5	LIST plist,*add_ptr;
176	gwlarson	3.1	{
177			int t_id,e_id0,e_id1,e_id2;
178			int v_id0,v_id1,v_id2;
179			LIST *lptr;
180
181			lptr = plist;
182			e_id0 = (int)LIST_DATA(lptr);
183			v_id0 = E_NTH_VERT(e_id0,0);
184			lptr = LIST_NEXT(lptr);
185			while(LIST_NEXT(lptr) != plist)
186			{
187			e_id1 = (int)LIST_DATA(lptr);
188			v_id1 = E_NTH_VERT(e_id1,0);
189			v_id2 = E_NTH_VERT(e_id1,1);
190			lptr = LIST_NEXT(lptr);
191
192	gwlarson	3.7	if(LIST_NEXT(lptr) != plist)
193			e_id2 = 0;
194	gwlarson	3.1	else
195			e_id2 = (int)LIST_DATA(lptr);
196	gwlarson	3.7	t_id = smTriangulate_add_tri(sm,v_id0,v_id1,v_id2,e_id0,e_id1,&e_id2);
197			add_ptr = push_data(add_ptr,t_id);
198	gwlarson	3.1	e_id0 = -e_id2;
199			}
200			free_list(plist);
201	gwlarson	3.3	return(TRUE);
202	gwlarson	3.1	}
203	gwlarson	3.7	#ifdef TEST_DRIVER
204			FVECT Norm[500],B_V[500];
205			int Ncnt,Bcnt,Del=0;
206	gwlarson	3.1	#endif
207
208	gwlarson	3.6
209	gwlarson	3.7	/* Triangulate the polygon defined by plist, and generating vertex p_id.
210			Return list of added triangles in list add_ptr. Returns TRUE if
211			successful, FALSE otherwise. This is NOT a general triangulation routine,
212			assumes polygon star relative to id
213			*/
214	gwlarson	3.6
215			int
216	gwlarson	3.7	smTriangulate(sm,id,plist,add_ptr)
217	gwlarson	3.6	SM *sm;
218			int id;
219	gwlarson	3.5	LIST plist,*add_ptr;
220	gwlarson	3.1	{
221	gwlarson	3.6	LIST l,prev,*t;
222			FVECT v0,v1,v2,n,p;
223	gwlarson	3.7	int is_tri,is_convex,cut,t_id,id0,id1,id2,e2,e1,enew;
224	gwlarson	3.6	double dp;
225	gwlarson	3.8	static int debug=0;
226	gwlarson	3.6
227	gwlarson	3.7	VSUB(p,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
228			enew = 0;
229			is_convex = TRUE;
230			cut = is_tri= FALSE;
231	gwlarson	3.6	l = prev = plist;
232	gwlarson	3.7
233			/* get v0,v1 */
234			e1 = (int)LIST_DATA(l);
235			id0 = E_NTH_VERT(e1,0);
236			id1 = E_NTH_VERT(e1,1);
237			VSUB(v0,SM_NTH_WV(sm,id0),SM_VIEW_CENTER(sm));
238			VSUB(v1,SM_NTH_WV(sm,id1),SM_VIEW_CENTER(sm));
239			#ifdef TEST_DRIVER
240			Del = TRUE;
241			VCOPY(B_V[0],v0);
242			VCOPY(B_V[1],v1);
243			Bcnt = 2;
244			Ncnt = 0;
245			#endif
246	gwlarson	3.6	while(l)
247			{
248			l = LIST_NEXT(l);
249	gwlarson	3.7	/* Get v2 */
250	gwlarson	3.6	e2 = (int)LIST_DATA(l);
251			id2 = E_NTH_VERT(e2,1);
252	gwlarson	3.7	VSUB(v2,SM_NTH_WV(sm,id2),SM_VIEW_CENTER(sm));
253			#ifdef TEST_DRIVER
254			VCOPY(B_V[Bcnt++],v2);
255			#endif
256	gwlarson	3.6	if(LIST_NEXT(LIST_NEXT(l)) == prev)
257	gwlarson	3.7	{/* Check if have a triangle */
258			is_tri = TRUE;
259			break;
260			}
261
262			/* determine if v0-v1-v2 is convex:defined clockwise on the sphere-
263			so switch orientation
264			*/
265			if(convex_angle(v2,v1,v0))
266	gwlarson	3.6	{
267	gwlarson	3.7	/* test if safe to cut off v0-v1-v2 by testing if p lies outside of
268			triangle v0-v1-v2: if so, because plist is the star polygon around p,
269			the new edge v2-v0 cannot intersect any existing edges
270			*/
271			VCROSS(n,v0,v2);
272			dp = DOT(n,p);
273			if(dp <= 0.0)
274			{
275			/* remove edges e1,e2 and add triangle id0,id1,id2 */
276			enew = 0;
277			t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew);
278			cut = TRUE;
279			add_ptr = push_data(add_ptr,t_id);
280			/* Insert edge enew into the list, reuse prev list element */
281			LIST_NEXT(prev) = LIST_NEXT(l);
282			LIST_DATA(prev) = e1 = -enew;
283			/* If removing head of list- reset plist pointer */
284			if(l== plist)
285			plist = prev;
286			/* free list element for e2 */
287			LIST_NEXT(l)=NULL;
288			free_list(l);
289			l = prev;
290			VCOPY(v1,v2);
291			id1 = id2;
292			continue;
293			}
294	gwlarson	3.6	}
295	gwlarson	3.7	else
296			is_convex = FALSE;
297			VCOPY(v0,v1);
298			VCOPY(v1,v2);
299			id0 = id1;
300			id1 = id2;
301			e1 = e2;
302			/* check if gone around circular list without adding any
303			triangles: prevent infinite loop */
304			if(l == plist)
305	gwlarson	3.6	{
306	gwlarson	3.7	if(LIST_NEXT(LIST_NEXT(l)) == prev)
307			{/* Check if have a triangle */
308			is_tri = TRUE;
309			break;
310			}
311
312			if(is_convex)
313	gwlarson	3.6	break;
314	gwlarson	3.7	if(!cut)
315			{
316			#ifdef DEBUG
317			eputs("smTriangulate():Unable to triangulate\n");
318			#endif
319	gwlarson	3.8	free_list(l);
320	gwlarson	3.7	while(*add_ptr)
321			{
322			t_id = pop_list(add_ptr);
323			smDelete_tri(sm,t_id);
324			}
325			return(FALSE);
326	gwlarson	3.8	}
327
328	gwlarson	3.7	cut = FALSE;
329			is_convex = TRUE;
330	gwlarson	3.6	}
331	gwlarson	3.7	prev = l;
332	gwlarson	3.6	}
333			if(is_tri)
334			{
335			l = LIST_NEXT(l);
336	gwlarson	3.7	enew = (int)LIST_DATA(l);
337			t_id = smTriangulate_add_tri(sm,id0,id1,id2,e1,e2,&enew);
338	gwlarson	3.6	add_ptr = push_data(add_ptr,t_id);
339			free_list(l);
340	gwlarson	3.7	}
341	gwlarson	3.6	else
342	gwlarson	3.7	if(!smTriangulateConvex(sm,l,add_ptr))
343	gwlarson	3.6	return(FALSE);
344
345	gwlarson	3.7	/* Set triangle adjacencies based on edge adjacencies */
346			FOR_ALL_EDGES(enew)
347	gwlarson	3.1	{
348	gwlarson	3.7	id0 = E_NTH_TRI(enew,0);
349			id1 = E_NTH_TRI(enew,1);
350	gwlarson	3.1
351	gwlarson	3.7	e1 = (T_WHICH_V(SM_NTH_TRI(sm,id0),E_NTH_VERT(enew,0))+2)%3;
352			T_NTH_NBR(SM_NTH_TRI(sm,id0),e1) = id1;
353
354			e2 = (T_WHICH_V(SM_NTH_TRI(sm,id1),E_NTH_VERT(enew,1))+2)%3;
355			T_NTH_NBR(SM_NTH_TRI(sm,id1),e2) = id0;
356	gwlarson	3.1	}
357	gwlarson	3.7	return(TRUE);
358	gwlarson	3.1	}
359
360			eIn_tri(e,t)
361			int e;
362			TRI *t;
363			{
364
365			if(T_NTH_V(t,0)==E_NTH_VERT(e,0))
366			return(T_NTH_V(t,1)==E_NTH_VERT(e,1)\|\|T_NTH_V(t,2)==E_NTH_VERT(e,1));
367			else
368			if(T_NTH_V(t,1)==E_NTH_VERT(e,0))
369			return(T_NTH_V(t,0)==E_NTH_VERT(e,1)\|\|T_NTH_V(t,2)==E_NTH_VERT(e,1));
370			else if(T_NTH_V(t,2)==E_NTH_VERT(e,0))
371			return(T_NTH_V(t,0)==E_NTH_VERT(e,1)\|\|T_NTH_V(t,1)==E_NTH_VERT(e,1));
372	gwlarson	3.7
373	gwlarson	3.1	return(FALSE);
374			}
375	gwlarson	3.6
376	gwlarson	3.7	/* Test the new set of triangles for Delaunay condition. 'Edges' contains
377			all of the new edges added. The CCW triangle assoc with each edge is
378			tested against the opposite vertex of the CW triangle. If the vertex
379			lies inside the circle defined by the CCW triangle- the edge is swapped
380			for the opposite diagonal
381			*/
382	gwlarson	3.8	smFixEdges(sm,add_list)
383	gwlarson	3.1	SM *sm;
384	gwlarson	3.5	LIST *add_list;
385	gwlarson	3.1	{
386	gwlarson	3.6	int e,t0_id,t1_id,e_new,e0,e1,e0_next,e1_next;
387			int i,v0_id,v1_id,v2_id,p_id,t0_nid,t1_nid;
388	gwlarson	3.1	FVECT v0,v1,v2,p,np,v;
389	gwlarson	3.7	TRI t0,t1;
390	gwlarson	3.5
391	gwlarson	3.1	FOR_ALL_EDGES(e)
392			{
393	gwlarson	3.6	t0_id = E_NTH_TRI(e,0);
394			t1_id = E_NTH_TRI(e,1);
395			if((t0_id==INVALID) \|\| (t1_id==INVALID))
396	gwlarson	3.1	{
397			#ifdef DEBUG
398	gwlarson	3.7	error(CONSISTENCY,"smFix_edges: Unassigned edge nbr\n");
399	gwlarson	3.1	#endif
400			}
401	gwlarson	3.7	t0 = SM_NTH_TRI(sm,t0_id);
402			t1 = SM_NTH_TRI(sm,t1_id);
403			e0 = T_NTH_NBR_PTR(t1_id,t0);
404			e1 = T_NTH_NBR_PTR(t0_id,t1);
405
406	gwlarson	3.6	v0_id = E_NTH_VERT(e,0);
407			v1_id = E_NTH_VERT(e,1);
408	gwlarson	3.7	v2_id = T_NTH_V(t0,e0);
409			p_id = T_NTH_V(t1,e1);
410	gwlarson	3.1
411	gwlarson	3.11	smDir(sm,v0,v0_id);
412			smDir(sm,v1,v1_id);
413			smDir(sm,v2,v2_id);
414	gwlarson	3.1
415	gwlarson	3.6	VCOPY(p,SM_NTH_WV(sm,p_id));
416	gwlarson	3.1	VSUB(p,p,SM_VIEW_CENTER(sm));
417			if(point_in_cone(p,v0,v1,v2))
418			{
419	gwlarson	3.8	smTris_swap_edge(sm,t0_id,t1_id,e0,e1,&t0_nid,&t1_nid,&add_list);
420	gwlarson	3.1
421	gwlarson	3.7	/* Adjust the triangle pointers of the remaining edges to be
422			processed
423			*/
424	gwlarson	3.1	FOR_ALL_EDGES_FROM(e,i)
425			{
426	gwlarson	3.6	if(E_NTH_TRI(i,0)==t0_id \|\| E_NTH_TRI(i,0)==t1_id)
427	gwlarson	3.1	{
428	gwlarson	3.6	if(eIn_tri(i,SM_NTH_TRI(sm,t0_nid)))
429			SET_E_NTH_TRI(i,0,t0_nid);
430	gwlarson	3.1	else
431	gwlarson	3.6	SET_E_NTH_TRI(i,0,t1_nid);
432	gwlarson	3.1	}
433
434	gwlarson	3.6	if(E_NTH_TRI(i,1)==t0_id \|\| E_NTH_TRI(i,1)==t1_id)
435	gwlarson	3.1	{
436	gwlarson	3.6	if(eIn_tri(i,SM_NTH_TRI(sm,t0_nid)))
437			SET_E_NTH_TRI(i,1,t0_nid);
438	gwlarson	3.1	else
439	gwlarson	3.6	SET_E_NTH_TRI(i,1,t1_nid);
440	gwlarson	3.1	}
441			}
442	gwlarson	3.6	t0_id = t0_nid;
443			t1_id = t1_nid;
444	gwlarson	3.1	e_new = eNew_edge();
445	gwlarson	3.6	SET_E_NTH_VERT(e_new,0,p_id);
446			SET_E_NTH_VERT(e_new,1,v2_id);
447			SET_E_NTH_TRI(e_new,0,t0_id);
448			SET_E_NTH_TRI(e_new,1,t1_id);
449	gwlarson	3.1	}
450			}
451	gwlarson	3.7	/* Add/Delete the appropriate triangles from the stree */
452	gwlarson	3.8	smUpdate_locator(sm,add_list);
453	gwlarson	3.1	}
454
455	gwlarson	3.7	/* Remove vertex "id" from the mesh- and retriangulate the resulting
456			hole: Returns TRUE if successful, FALSE otherwise.
457			*/
458	gwlarson	3.1	int
459	gwlarson	3.7	smRemoveVertex(sm,id)
460	gwlarson	3.1	SM *sm;
461			int id;
462			{
463	gwlarson	3.8	LIST b_list,add_list,*del_list;
464			int t_id,i;
465			static int cnt=0;
466			OBJECT optr,os;
467	gwlarson	3.7	/* generate list of edges that form the boundary of the
468			polygon formed by the triangles adjacent to vertex 'id'
469	gwlarson	3.1	*/
470	gwlarson	3.8	del_list = NULL;
471	gwlarson	3.11	b_list = smVertexStar(sm,id,&del_list);
472	gwlarson	3.10	if(!b_list)
473			{
474			if(del_list)
475			free_list(del_list);
476			return(FALSE);
477			}
478	gwlarson	3.5	add_list = NULL;
479	gwlarson	3.7	/* Triangulate polygonal hole */
480			if(!smTriangulate(sm,id,b_list,&add_list))
481	gwlarson	3.6	{
482	gwlarson	3.8	free_list(del_list);
483	gwlarson	3.6	return(FALSE);
484			}
485	gwlarson	3.8	else
486			{
487			while(del_list)
488			{
489			t_id = pop_list(&del_list);
490			smDelete_tri(sm,t_id);
491			}
492			}
493	gwlarson	3.7	/* Fix up new triangles to be Delaunay-delnode contains set of
494			triangles to delete,add_list is the set of new triangles to add
495			*/
496	gwlarson	3.8	smFixEdges(sm,add_list);
497	gwlarson	3.7
498	gwlarson	3.1	return(TRUE);
499			}
500
501
502
503	gwlarson	3.8
504	gwlarson	3.1
505
506
507
508
509
510
511
512
513
514