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;

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

  tri = SM_NTH_TRI(sm,id);
  /* Add to the free_list */
  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
  */
  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);

  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)
        error(CONSISTENCY,"Too many edges in vertex loop\n");
      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 
*smVertexPolygon(sm,id,delptr)
SM *sm;
int id;
QUADTREE *delptr;
{
    TRI *tri,*t_next;
    LIST *elist,*end;
    int e,t_id,v_next,t_next_id,b_id,v_id;
    OBJECT del_set[2];

    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;

    /* 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,b_id)) != t_id)
    {
      e = eNew_edge();
      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);
      qtaddelem(*delptr,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;

    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,delptr)
   SM *sm;
   LIST *add_list;
   QUADTREE *delptr;

{
    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_in_cone(sm,v0,v0_id);
        smDir_in_cone(sm,v1,v1_id);
        smDir_in_cone(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,
                            delptr);
            
           /* 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,qtqueryset(*delptr));

}

/* 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;
    QUADTREE delnode=-1;
    int t_id;

    /* generate list of edges that form the boundary of the
       polygon formed by the triangles adjacent to vertex 'id'
     */
    b_list = smVertexPolygon(sm,id,&delnode);

    add_list = NULL;
    /* Triangulate polygonal hole  */
    if(!smTriangulate(sm,id,b_list,&add_list))
    {
      qtfreeleaf(delnode);
      return(FALSE);
    }
    /* 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,&delnode);


    qtfreeleaf(delnode);
    return(TRUE);
}
   

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