ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/sm.c
(Generate patch)

Comparing ray/src/hd/sm.c (file contents):
Revision 3.14 by gwlarson, Fri Mar 5 16:33:17 1999 UTC vs.
Revision 3.18 by greg, Fri Jun 20 00:25:49 2003 UTC

# Line 1 | Line 1
1 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ SGI";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6
4   /*
5   *  sm.c
6   */
# Line 11 | Line 8 | static char SCCSid[] = "$SunId$ SGI";
8   #include "sm_flag.h"
9   #include "sm_list.h"
10   #include "sm_geom.h"
14 #include "sm_qtree.h"
15 #include "sm_stree.h"
11   #include "sm.h"
12  
18
13   SM *smMesh = NULL;
14   double smDist_sum=0;
21 #define S_INC 1000
22 int smNew_tri_cnt=0;
23 int smNew_tri_size=0;
15  
16 < #define SM_MIN_SAMP_DIFF 1e-4/* min edge length in radians */
17 <
18 < /* Each edge extends .5536 radians - 31.72 degrees */
16 > #ifdef DEBUG
17 > int Tcnt=0,Wcnt=0;
18 > #endif
19 > /* Each edge extends .372 radians  */
20   static FVECT icosa_verts[162] =
21   {{-0.018096,0.495400,0.868477},{0.018614,-0.554780,0.831789},
22   {0.850436,0.011329,0.525956},{0.850116,0.048016,-0.524402},
# Line 131 | Line 123 | static int icosa_tris[320][3] =
123   {34,112,114},{112,33,113},{114,113,1},{112,113,114},{32,111,107},{111,33,112},
124   {107,112,34},{111,112,107},{9,109,116},{109,32,115},{116,115,36},{109,115,116},
125   {32,106,118},{106,8,117},{118,117,35},{106,117,118},{36,119,121},{119,35,120},
126 < {121,120,11},{119,120,121},{32,118,115},{118,35,119},{115,119,36},{118,119,115},
127 < {10,122,124},{122,37,123},{124,123,39},{122,123,124},{37,125,127},{125,11,126},
128 < {127,126,38},{125,126,127},{39,128,130},{128,38,129},{130,129,3},{128,129,130},
126 > {121,120,11},{119,120,121},{32,118,115},{118,35,119},{115,119,36},
127 > {118,119,115},{10,122,124},{122,37,123},{124,123,39},{122,123,124},
128 > {37,125,127},{125,11,126},{127,126,38},{125,126,127},{39,128,130},
129 > {128,38,129},{130,129,3},{128,129,130},
130   {37,127,123},{127,38,128},{123,128,39},{127,128,123},{11,125,132},{125,37,131},
131 < {132,131,41},{125,131,132},{37,122,134},{122,10,133},{134,133,40},{122,133,134},
132 < {41,135,137},{135,40,136},{137,136,5},{135,136,137},{37,134,131},{134,40,135},
131 > {132,131,41},{125,131,132},{37,122,134},{122,10,133},{134,133,40},
132 > {122,133,134},{41,135,137},{135,40,136},{137,136,5},{135,136,137},
133 > {37,134,131},{134,40,135},
134   {131,135,41},{134,135,131},{6,65,94},{65,18,138},{94,138,28},{65,138,94},
135 < {18,62,139},{62,2,49},{139,49,13},{62,49,139},{28,140,97},{140,13,46},{97,46,0},
136 < {140,46,97},{18,139,138},{139,13,140},{138,140,28},{139,140,138},{1,44,114},
135 > {18,62,139},{62,2,49},{139,49,13},{62,49,139},{28,140,97},{140,13,46},
136 > {97,46,0},{140,46,97},{18,139,138},{139,13,140},{138,140,28},{139,140,138},
137 > {1,44,114},
138   {44,14,141},{114,141,34},{44,141,114},{14,50,142},{50,2,68},{142,68,21},
139   {50,68,142},{34,143,108},{143,21,73},{108,73,8},{143,73,108},{14,142,141},
140   {142,21,143},{141,143,34},{142,143,141},{0,52,98},{52,16,144},{98,144,29},
# Line 155 | Line 150 | static int icosa_tris[320][3] =
150   {155,38,126},{120,126,11},{155,126,120},{20,154,153},{154,38,155},{153,155,35},
151   {154,155,153},{7,81,101},{81,23,156},{101,156,30},{81,156,101},{23,78,157},
152   {78,5,136},{157,136,40},{78,136,157},{30,158,104},{158,40,133},{104,133,10},
153 < {158,133,104},{23,157,156},{157,40,158},{156,158,30},{157,158,156},{11,132,121},
153 > {158,133,104},{23,157,156},{157,40,158},{156,158,30},{157,158,156},{
154 > 11,132,121},
155   {132,41,159},{121,159,36},{132,159,121},{41,137,160},{137,5,84},{160,84,26},
156   {137,84,160},{36,161,116},{161,26,89},{116,89,9},{161,89,116},{41,160,159},
157   {160,26,161},{159,161,36},{160,161,159}};
# Line 195 | Line 191 | static int icosa_nbrs[320][3] =
191   {182,180,181},{187,308,190},{294,187,189},{293,309,187},{186,184,185},
192   {191,177,180},{185,191,182},{184,178,191},{190,188,189},{195,101,42},
193   {204,195,41},{206,102,195},{194,192,193},{199,204,38},{10,199,37},{9,205,199},
194 < {198,196,197},{203,105,206},{6,203,205},{5,106,203},{202,200,201},{207,193,196},
195 < {201,207,198},{200,194,207},{206,204,205},{211,138,0},{220,211,2},{222,136,211},
196 < {210,208,209},{215,220,8},{48,215,10},{50,221,215},{214,212,213},{219,130,222},
194 > {198,196,197},{203,105,206},{6,203,205},{5,106,203},{202,200,201},
195 > {207,193,196},{201,207,198},{200,194,207},{206,204,205},{211,138,0},
196 > {220,211,2}, {222,136,211},{210,208,209},{215,220,8},{48,215,10},{50,221,215},
197 > {214,212,213},{219,130,222},
198   {56,219,221},{58,128,219},{218,216,217},{223,209,212},{217,223,214},
199   {216,210,223},{222,220,221},{227,106,16},{236,227,18},{238,104,227},
200 < {226,224,225},{231,236,24},{64,231,26},{66,237,231},{230,228,229},{235,98,238},
200 > {226,224,225},{231,236,24},{64,231,26},{66,237,231},{230,228,229},
201 > {235,98,238},
202   {72,235,237},{74,96,235},{234,232,233},{239,225,228},{233,239,230},
203   {232,226,239},{238,236,237},{243,133,90},{252,243,89},{254,134,243},
204 < {242,240,241},{247,252,86},{26,247,85},{25,253,247},{246,244,245},{251,137,254},
204 > {242,240,241},{247,252,86},{26,247,85},{25,253,247},{246,244,245},
205 > {251,137,254},
206   {22,251,253},{21,138,251},{250,248,249},{255,241,244},{249,255,246},
207   {248,242,255},{254,252,253},{259,122,160},{268,259,162},{270,120,259},
208   {258,256,257},{263,268,168},{32,263,170},{34,269,263},{262,260,261},
# Line 233 | Line 232 | FVECT FrustumNear[4],FrustumFar[4];
232   double dev_zmin=.01,dev_zmax=1000;
233   #endif
234  
235 +
236 + char *
237 + tempbuf(len,resize)                     /* get a temporary buffer */
238 + unsigned  len;
239 + int resize;
240 + {
241 +  static char  *tempbuf = NULL;
242 +  static unsigned  tempbuflen = 0;
243 +
244 + #ifdef DEBUG
245 +        static int in_use=FALSE;
246 +
247 +        if(len == -1)
248 +          {
249 +            in_use = FALSE;
250 +            return(NULL);
251 +          }
252 +        if(!resize && in_use)
253 +        {
254 +            eputs("Buffer in use:cannot allocate:tempbuf()\n");
255 +            return(NULL);
256 +        }
257 + #endif
258 +        if (len > tempbuflen) {
259 +                if (tempbuflen > 0)
260 +                        tempbuf = realloc((void *)tempbuf, len);
261 +                else
262 +                        tempbuf = malloc(len);
263 +                tempbuflen = tempbuf==NULL ? 0 : len;
264 +        }
265 + #ifdef DEBUG
266 +        in_use = TRUE;
267 + #endif
268 +        return(tempbuf);
269 + }
270 +
271   smDir(sm,ps,id)
272     SM *sm;
273     FVECT ps;
274 <   int id;
274 >   S_ID id;
275   {
276      VSUB(ps,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
277      normalize(ps);
# Line 261 | Line 296 | smInit_mesh(sm)
296   {
297      /* Reset the triangle counters */
298      SM_NUM_TRI(sm) = 0;
264    SM_SAMPLE_TRIS(sm) = 0;
299      SM_FREE_TRIS(sm) = -1;
266    SM_AVAILABLE_TRIS(sm) = -1;
300      smClear_flags(sm,-1);
301   }
302  
# Line 274 | Line 307 | smInit_mesh(sm)
307   smClear(sm)
308   SM *sm;
309   {
310 <  smClear_samples(sm);
311 <  smClear_locator(sm);
310 >  smDist_sum = 0;
311 >  smInit_samples(sm);
312 >  smInit_locator(sm);
313    smInit_mesh(sm);
314   }
315  
282 static int realloc_cnt=0;
283
284 int
285 smRealloc_mesh(sm)
286 SM *sm;
287 {
288  int fbytes,i,max_tris,max_verts;
289  
290  if(realloc_cnt <=0)
291    realloc_cnt = SM_MAX_TRIS(sm);
292
293  max_tris = SM_MAX_TRIS(sm) + realloc_cnt;
294  fbytes = FLAG_BYTES(max_tris);
295
296  if(!(SM_TRIS(sm) = (TRI *)realloc(SM_TRIS(sm),max_tris*sizeof(TRI))))
297    goto memerr;
298
299  for(i=0; i< T_FLAGS; i++)
300   if(!(SM_NTH_FLAGS(sm,i)=(int4 *)realloc((char *)SM_NTH_FLAGS(sm,i),fbytes)))
301        goto memerr;
302
303  SM_MAX_TRIS(sm) = max_tris;
304  return(max_tris);
305
306 memerr:
307        error(SYSTEM, "out of memory in smRealloc_mesh()");
308 }
309
316   /* Allocate and initialize a new mesh with max_verts and max_tris */
317   int
318   smAlloc_mesh(sm,max_verts,max_tris)
# Line 317 | Line 323 | int max_verts,max_tris;
323  
324      fbytes = FLAG_BYTES(max_tris);
325  
326 <    if(!(SM_TRIS(sm) = (TRI *)realloc(NULL,max_tris*sizeof(TRI))))
326 >    if(!(SM_TRIS(sm) = (TRI *)malloc(max_tris*sizeof(TRI))))
327        goto memerr;
328  
323    if(!(SM_VERTS(sm) = (VERT *)realloc(NULL,max_verts*sizeof(VERT))))
324      goto memerr;
325
329      for(i=0; i< T_FLAGS; i++)
330 <      if(!(SM_NTH_FLAGS(sm,i)=(int4 *)realloc(NULL,fbytes)))
330 >      if(!(SM_NTH_FLAGS(sm,i)=(int32 *)malloc(fbytes)))
331          goto memerr;
332  
333      SM_MAX_VERTS(sm) = max_verts;
334      SM_MAX_TRIS(sm) = max_tris;
335  
333    realloc_cnt = max_verts >> 1;
334
336      smInit_mesh(sm);
337  
338      return(max_tris);
# Line 340 | Line 341 | memerr:
341   }
342  
343  
344 < int
345 < compress_set(os)
346 < OBJECT *os;
347 < {
348 <  int i,j;
349 <
349 <  for (i = 1, j = os[0]; i <= j; i++)
350 <  {
351 <    if(SM_T_ID_VALID(smMesh,os[i]))
352 <      continue;
353 <    if(i==j)
354 <      break;
355 <    while(!SM_T_ID_VALID(smMesh,os[j]) && (i < j))
356 <            j--;
357 <    if(i==j)
358 <      break;
359 <    os[i] = os[j--];
360 <  }
361 <  i--;
362 <        
363 <  os[0] = i;
364 <  return(i);
365 <
366 < }
367 <
344 > /*
345 > * smAlloc_tri(sm) : Allocate a new mesh triangle
346 > * SM *sm;                    : mesh
347 > *
348 > *  Returns ptr to next available tri
349 > */
350   int
351   smAlloc_tri(sm)
352   SM *sm;
# Line 375 | Line 357 | SM *sm;
357    /* Otherwise, have we reached the end of the list yet?*/
358    if(SM_NUM_TRI(sm) < SM_MAX_TRIS(sm))
359      return(SM_NUM_TRI(sm)++);
378  if((id = SM_AVAILABLE_TRIS(sm)) != -1)
379  {
380    SM_AVAILABLE_TRIS(sm) =  T_NEXT_AVAILABLE(SM_NTH_TRI(sm,id));
381    return(id);
382  }
360    if((id = SM_FREE_TRIS(sm)) != -1)
361    {
362 <    dosets(compress_set);
386 <    SM_AVAILABLE_TRIS(sm) = T_NEXT_FREE(SM_NTH_TRI(sm,id));
387 <    SM_FREE_TRIS(sm) = -1;
362 >    SM_FREE_TRIS(sm) =  T_NEXT_FREE(SM_NTH_TRI(sm,id));
363      return(id);
364    }
365 <  /*Else: realloc */
366 <  smRealloc_mesh(sm);
367 <  return(SM_NUM_TRI(sm)++);
365 >
366 >  error(CONSISTENCY,"smAlloc_tri: no more available triangles\n");
367 >  return(INVALID);
368   }
369  
370 +
371   smFree_mesh(sm)
372   SM *sm;
373   {
# Line 399 | Line 375 | SM *sm;
375  
376    if(SM_TRIS(sm))
377      free(SM_TRIS(sm));
402  if(SM_VERTS(sm))
403    free(SM_VERTS(sm));
378    for(i=0; i< T_FLAGS; i++)
379      if(SM_NTH_FLAGS(sm,i))
380        free(SM_NTH_FLAGS(sm,i));
381   }
382  
383 <  
384 < /* Initialize/clear global smL sample list for at least n samples */
383 >
384 > /*
385 > * smAlloc(max_samples) : Initialize/clear global sample list for at least
386 > *                       max_samples
387 > * int max_samples;    
388 > *
389 > */
390   smAlloc(max_samples)
391     register int max_samples;
392   {
# Line 434 | Line 413 | smAlloc(max_samples)
413    /* First allocate at least n samples + extra points:at least enough
414     necessary to form the BASE MESH- Default = 4;
415    */
416 <  SM_SAMP(smMesh) = sAlloc(&n,SM_EXTRA_POINTS);
416 >  SM_SAMP(smMesh) = sAlloc(&n,SM_BASE_POINTS);
417  
418 <  total_points = n + SM_EXTRA_POINTS;
418 >  total_points = n + SM_BASE_POINTS;
419  
420 <  max_tris = total_points*2;
420 >  /* Need enough tris for base mesh, + 2 for each sample, plus extra
421 >     since triangles are not deleted until added and edge flips (max 4)
422 >  */
423 >  max_tris = n*2 + SM_BASE_TRIS + 10;
424    /* Now allocate space for mesh vertices and triangles */
425    max_tris = smAlloc_mesh(smMesh, total_points, max_tris);
426  
# Line 453 | Line 435 | smInit_sm(sm,vp)
435   SM *sm;
436   FVECT vp;
437   {
456
457  smDist_sum = 0;
458  smNew_tri_cnt = 0;
459
438    VCOPY(SM_VIEW_CENTER(sm),vp);
439 <  smInit_locator(sm,vp);
462 <  smInit_samples(sm);
463 <  smInit_mesh(sm);  
439 >  smClear(sm);
440    smCreate_base_mesh(sm,SM_DEFAULT);
441   }
442  
# Line 474 | Line 450 | FVECT vp;
450   * If n is <= 0, then clear data structures.  Returns number samples
451   * actually allocated.
452   */
477
453   int
454   smInit(n)
455     register int n;
456   {
457    int max_vertices;
458 <
458 >  sleep(5);
459    /* If n <=0, Just clear the existing structures */
460    if(n <= 0)
461    {
462 + #ifdef DEBUG
463 +  fprintf(stderr,"Tcnt=%d Wcnt = %d, avg = %f\n",Tcnt,Wcnt,(float)Wcnt/Tcnt);
464 + #endif
465      smClear(smMesh);
466      return(0);
467    }
# Line 491 | Line 469 | smInit(n)
469    /* Total mesh vertices includes the sample points and the extra vertices
470       to form the base mesh
471    */
472 <  max_vertices = n + SM_EXTRA_POINTS;
472 >  max_vertices = n + SM_BASE_POINTS;
473  
474    /* If the current mesh contains enough room, clear and return */
475    if(smMesh && (max_vertices <= SM_MAX_VERTS(smMesh)) && SM_MAX_SAMP(smMesh) <=
# Line 505 | Line 483 | smInit(n)
483    */
484    smAlloc(n);
485  
486 + #ifdef DEBUG
487 +  fprintf(stderr,"smInit: allocating room for %d samples\n",
488 +          SM_MAX_SAMP(smMesh));
489 + #endif
490    return(SM_MAX_SAMP(smMesh));
491   }
492  
# Line 527 | Line 509 | smLocator_apply(sm,v0,v1,v2,func,n)
509  
510   }
511  
512 + /*
513 + * QUADTREE
514 + * insert_samp(argptr,root,qt,parent,q0,q1,q2,bc,scale,rev,f,n,doneptr)
515 + *     Callback function called from quadtree traversal when leaf is reached
516 + *   int *argptr;                 :ptr to sample id number
517 + *   int root;                    :quadtree root from which traversal started
518 + *   QUADTREE qt,parent;          :current quadtree node and its parent
519 + *   BCOORD q0[3],q1[3],q2[3],bc[3]; :barycentric coordinates of the current
520 + *                                 quadtree node (q0,q1,q2) and sample (bc)
521 + *   unsigned int scale,rev;      :scale factor relative to which level at in
522 + *                                 tree, rev indicates if traversed of child 3
523 + *   FUNC f;                      :function structure so can recurse
524 + *   int n,*doneptr;              :n indicates which level at in quadtree, and
525 + *                                 doneptr can be set to indicate whether nbr
526 + *                                 sample has been found
527 + */
528   QUADTREE
529 < insert_tri(argptr,root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
532 <               s0,s1,s2,sq0,sq1,sq2,rev,f,n)
529 > insert_samp(argptr,root,qt,parent,q0,q1,q2,bc,scale,rev,f,n,doneptr)
530       int *argptr;
531       int root;
532 <     QUADTREE qt;
533 <     BCOORD q0[3],q1[3],q2[3];
534 <     BCOORD t0[3],t1[3],t2[3];
538 <     BCOORD dt10[3],dt21[3],dt02[3];
539 <     unsigned int scale,s0,s1,s2,sq0,sq1,sq2,rev;
532 >     QUADTREE qt,parent;
533 >     BCOORD q0[3],q1[3],q2[3],bc[3];
534 >     unsigned int scale,rev;
535       FUNC f;
536 <     int n;
536 >     int n,*doneptr;
537   {
538 <  OBJECT *optr,*tptr,t_set[QT_MAXSET+1];
539 <  int i,t_id;
545 <  TRI *t;
546 <  FVECT tri[3];
547 <  BCOORD b0[3],b1[3],b2[3];
548 <  BCOORD tb0[3],tb1[3],tb2[3];
549 <  BCOORD db10[3],db21[3],db02[3];
550 <  BCOORD a,b,c,qa,qb,qc;
551 <  FUNC tfunc;
538 >  S_ID s_id;
539 >  S_ID *optr;
540  
541 <  t_id = *argptr;
554 <
541 >  s_id = ((S_ARGS *)argptr)->s_id;
542    /* If the set is empty - just add */
543    if(QT_IS_EMPTY(qt))
544 <    return(qtadduelem(qt,t_id));
545 <
544 >  {
545 >    qt = qtadduelem(qt,s_id);
546 >    SM_S_NTH_QT(smMesh,s_id) = qt;
547 >    if(((S_ARGS *)argptr)->n_id)
548 >      *doneptr = FALSE;
549 >    else
550 >      *doneptr = TRUE;
551 >    return(qt);
552 >  }
553    /* If the set size is below expansion threshold,or at maximum
554       number of quadtree levels : just add */
555    optr = qtqueryset(qt);
556    if(QT_SET_CNT(optr) < QT_SET_THRESHOLD || (n > (QT_MAX_LEVELS-2)))
557 <    return(qtadduelem(qt,t_id));
558 <  /* otherwise: expand node- and insert tri, and reinsert existing tris
557 >  {
558 >    qt = qtadduelem(qt,s_id);
559 >    
560 >    SM_S_NTH_QT(smMesh,s_id) = qt;
561 >    if(((S_ARGS *)argptr)->n_id)
562 >      if(QT_SET_CNT(qtqueryset(qt)) > 1)
563 >      {
564 >        ((S_ARGS *)argptr)->n_id = qtqueryset(qt)[1];
565 >        *doneptr = TRUE;
566 >      }
567 >      else
568 >        *doneptr = FALSE;
569 >    else
570 >      *doneptr = TRUE;
571 >    return(qt);
572 >  }
573 >  /* otherwise: expand node- and insert sample, and reinsert existing samples
574       in set to children of new node
575    */
576 <  /* First tri compressing set */
577 < #ifdef NEWSETS0
578 <  if(qtcompressuelem(qt,compress_set) < QT_SET_THRESHOLD)
579 <    /* Get set: If greater than standard size: allocate memory to hold */
580 <    return(qtadduelem(qt,t_id));
581 < #endif
582 <    if(QT_SET_CNT(optr) > QT_MAXSET)
583 <    {
584 <       if(!(tptr = (OBJECT *)malloc((QT_SET_CNT(optr)+1)*sizeof(OBJECT))))
585 <        goto memerr;
586 <    }
587 <    else
588 <     tptr = t_set;
589 <    qtgetset(tptr,qt);
576 >  {
577 >      OBJECT *sptr,s_set[QT_MAXSET+1];
578 >      FUNC sfunc;
579 >      S_ARGS args;
580 >      FVECT p;
581 >      BCOORD bp[3];
582 >      int i;
583 >      
584 >      if(QT_SET_CNT(optr) > QT_MAXSET)
585 >      {
586 >          if(!(sptr = (OBJECT *)malloc((QT_SET_CNT(optr)+1)*sizeof(OBJECT))))
587 >             goto memerr;
588 >      }
589 >      else
590 >         sptr = s_set;
591  
592 <  /* subdivide node */
593 <  qtfreeleaf(qt);
594 <  qtSubdivide(qt);
592 >      qtgetset(sptr,qt);
593 >      
594 >      /* subdivide node */
595 >      qtfreeleaf(qt);
596 >      qtSubdivide(qt);
597 >      
598 >      /* Now add in all of the rest; */
599 >      F_FUNC(sfunc) = F_FUNC(f);
600 >      F_ARGS(sfunc) = (int *) (&args);
601 >      args.n_id = 0;
602 >      for(optr = sptr,i=QT_SET_CNT(sptr); i > 0; i--)
603 >      {
604 >          s_id = QT_SET_NEXT_ELEM(optr);
605 >          args.s_id = s_id;
606 >          VSUB(p,SM_NTH_WV(smMesh,s_id),SM_VIEW_CENTER(smMesh));
607 >          normalize(p);
608 >          vert_to_qt_frame(i,p,bp);    
609 >          if(rev)
610 >             qt= qtInsert_point_rev(root,qt,parent,q0,q1,q2,bp,
611 >                                    scale,sfunc,n,doneptr);
612 >          else
613 >             qt= qtInsert_point(root,qt,parent,q0,q1,q2,bp,
614 >                                scale,sfunc,n,doneptr);
615 >      }
616 >      /* Add current sample: have all of the information */
617 >      if(rev)
618 >         qt =qtInsert_point_rev(root,qt,parent,q0,q1,q2,bc,scale,f,n,doneptr);
619 >      else
620 >         qt = qtInsert_point(root,qt,parent,q0,q1,q2,bc,scale,f,n,doneptr);
621 >      
622 >      /* If we allocated memory: free it */
623  
624 <  a = q1[0];  b = q0[1]; c = q0[2];
625 <  qa = q0[0]; qb = q1[1];  qc = q2[2];
588 <  /* First add current tri: have all of the information */
589 <  if(rev)
590 <    qt = qtInsert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
591 <               s0,s1,s2,sq0,sq1,sq2,f,n);
592 <  else
593 <    qt = qtInsert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
594 <               s0,s1,s2,sq0,sq1,sq2,f,n);
624 >      if( sptr != s_set)
625 >         free(sptr);
626  
627 <  /* Now add in all of the rest; */
597 <  F_FUNC(tfunc) = F_FUNC(f);
598 <  for(optr = QT_SET_PTR(tptr),i=QT_SET_CNT(tptr); i > 0; i--)
599 <  {
600 <    t_id = QT_SET_NEXT_ELEM(optr);
601 <    t = SM_NTH_TRI(smMesh,t_id);
602 <    if(!T_IS_VALID(t))
603 <      continue;
604 <    F_ARGS(tfunc) = &t_id;
605 <    VSUB(tri[0],SM_T_NTH_WV(smMesh,t,0),SM_VIEW_CENTER(smMesh));
606 <    VSUB(tri[1],SM_T_NTH_WV(smMesh,t,1),SM_VIEW_CENTER(smMesh));
607 <    VSUB(tri[2],SM_T_NTH_WV(smMesh,t,2),SM_VIEW_CENTER(smMesh));
608 <    convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02);
609 <
610 <    /* Calculate appropriate sidedness values */
611 <    SIDES_GTR(b0,b1,b2,s0,s1,s2,a,b,c);
612 <    SIDES_GTR(b0,b1,b2,sq0,sq1,sq2,qa,qb,qc);
613 <    if(rev)
614 <      qt = qtInsert_tri_rev(root,qt,q0,q1,q2,b0,b1,b2,db10,db21,db02,scale,
615 <               s0,s1,s2,sq0,sq1,sq2,tfunc,n);
616 <    else
617 <      qt = qtInsert_tri(root,qt,q0,q1,q2,b0,b1,b2,db10,db21,db02,scale,
618 <               s0,s1,s2,sq0,sq1,sq2,tfunc,n);
627 >      return(qt);
628    }
620  /* If we allocated memory: free it */
621  if( tptr != t_set)
622    free(tptr);
623
624  return(qt);
629   memerr:
630      error(SYSTEM,"expand_node():Unable to allocate memory");
631  
632   }
633  
634  
635 <
636 < smLocator_add_tri(sm,t_id,v0_id,v1_id,v2_id)
637 < SM *sm;
638 < int t_id;
639 < int v0_id,v1_id,v2_id;
640 < {
641 <  FVECT tri[3];
642 <  FUNC f;
643 <
644 < #ifdef DEBUG
645 < if(T_NTH_NBR(SM_NTH_TRI(sm,t_id),0)== -1 ||
642 <   T_NTH_NBR(SM_NTH_TRI(sm,t_id),1)== -1 ||
643 <   T_NTH_NBR(SM_NTH_TRI(sm,t_id),2)== -1)
644 <  error("Invalid tri: smLocator_add_tri()\n");
645 <
646 < if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,t_id),0))) ||
647 <   !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,t_id),1))) ||
648 <   !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,t_id),2))))
649 <  error("Invalid tri: smLocator_add_tri()\n");
650 < #endif
651 <
652 <  VSUB(tri[0],SM_NTH_WV(sm,v0_id),SM_VIEW_CENTER(sm));
653 <  VSUB(tri[1],SM_NTH_WV(sm,v1_id),SM_VIEW_CENTER(sm));
654 <  VSUB(tri[2],SM_NTH_WV(sm,v2_id),SM_VIEW_CENTER(sm));
655 <
656 <  F_FUNC(f) = insert_tri;
657 <  F_ARGS(f) = &t_id;
658 <  stInsert_tri(SM_LOCATOR(sm),tri,f);
659 < }
660 <
661 < /* Add a triangle to the base array with vertices v1-v2-v3 */
635 > /*
636 > * int
637 > * smAdd_tri(sm,v0_id,v1_id,v2_id,tptr)
638 > *             : Add a triangle to the base array with vertices v0-v1-v2,
639 > *               returns ptr to triangle in tptr
640 > * SM *sm;                   : mesh
641 > * S_ID v0_id,v1_id,v2_id;    : sample ids of triangle vertices
642 > * TRI **tptr;               : ptr to set to triangle
643 > *
644 > *  Allocates and initializes mesh triangle with vertices specified.
645 > */
646   int
647 < smAdd_tri(sm, v0_id,v1_id,v2_id)
647 > smAdd_tri(sm, v0_id,v1_id,v2_id,tptr)
648   SM *sm;
649 < int v0_id,v1_id,v2_id;
649 > S_ID v0_id,v1_id,v2_id;
650 > TRI **tptr;
651   {
652      int t_id;
653      TRI *t;
654 +
655   #ifdef DEBUG
656      if(v0_id==v1_id || v0_id==v2_id || v1_id==v2_id)
657        error(CONSISTENCY,"smAdd_tri: invalid vertex ids\n");
672 #endif    
673    t_id = smAlloc_tri(sm);
658  
659 <    if(t_id == -1)
660 <      return(t_id);
659 > #endif
660 >    t_id = smAlloc_tri(sm);  
661 > #ifdef DEBUG
662 > #if DEBUG > 1
663 >    {
664 >      double v0[3],v1[3],v2[3],n[3];
665 >      double area,dp;
666 >      double ab[3],ac[3];
667 >      VSUB(v0,SM_NTH_WV(sm,v0_id),SM_VIEW_CENTER(sm));
668 >      VSUB(v1,SM_NTH_WV(sm,v1_id),SM_VIEW_CENTER(sm));
669 >      VSUB(v2,SM_NTH_WV(sm,v2_id),SM_VIEW_CENTER(sm));
670 >      normalize(v0);
671 >      normalize(v1);
672 >      normalize(v2);
673  
674 <    t = SM_NTH_TRI(sm,t_id);
674 >      VSUB(ab,v1,v2);
675 >      normalize(ab);
676 >      VSUB(ac,v0,v2);
677 >      normalize(ac);
678 >      VCROSS(n,ab,ac);
679 >      area = normalize(n);
680 >      if((dp=DOT(v0,n)) < 0.0)
681 >      {
682 >        eputs("backwards tri\n");
683 >      }
684 >    }
685 > #endif
686 > #endif
687  
688 +    t = SM_NTH_TRI(sm,t_id);
689 + #ifdef DEBUG
690      T_CLEAR_NBRS(t);
691 + #endif
692      /* set the triangle vertex ids */
693      T_NTH_V(t,0) = v0_id;
694      T_NTH_V(t,1) = v1_id;
# Line 689 | Line 700 | int v0_id,v1_id,v2_id;
700  
701      if(SM_BASE_ID(sm,v0_id) || SM_BASE_ID(sm,v1_id) || SM_BASE_ID(sm,v2_id))
702        SM_SET_NTH_T_BASE(sm,t_id);
703 <
703 >    else
704 >      SM_CLR_NTH_T_BASE(sm,t_id);
705      if(SM_DIR_ID(sm,v0_id) || SM_DIR_ID(sm,v1_id) || SM_DIR_ID(sm,v2_id))
706        SM_SET_NTH_T_BG(sm,t_id);
707 <
707 >    else
708 >      SM_CLR_NTH_T_BG(sm,t_id);
709      S_SET_FLAG(T_NTH_V(t,0));
710      S_SET_FLAG(T_NTH_V(t,1));
711      S_SET_FLAG(T_NTH_V(t,2));
# Line 700 | Line 713 | int v0_id,v1_id,v2_id;
713      SM_SET_NTH_T_ACTIVE(sm,t_id);
714      SM_SET_NTH_T_NEW(sm,t_id);
715  
716 <    SM_SAMPLE_TRIS(sm)++;
704 <    smNew_tri_cnt++;
705 <
716 >    *tptr = t;
717      /* return initialized triangle */
718      return(t_id);
719   }
720  
721  
722 < void
723 < smTris_swap_edge(sm,t_id,t1_id,e,e1,tn_id,tn1_id,add_ptr)
724 <   SM *sm;
725 <   int t_id,t1_id;
726 <   int e,e1;
727 <   int *tn_id,*tn1_id;
717 <   LIST **add_ptr;
722 >
723 > smRestore_Delaunay(sm,a,b,c,t,t_id,a_id,b_id,c_id)
724 > SM *sm;
725 > FVECT a,b,c;
726 > TRI *t;
727 > int t_id,a_id,b_id,c_id;
728   {
729 <    int verts[3],enext,eprev,e1next,e1prev;
730 <    TRI *n,*ta,*tb,*t,*t1;
731 <    FVECT p1,p2,p3;
732 <    int ta_id,tb_id;
723 <    /* form new diagonal (e relative to t, and e1 relative to t1)
724 <      defined by quadrilateral formed by t,t1- swap for the opposite diagonal
725 <   */
726 <    enext = (e+1)%3;
727 <    eprev = (e+2)%3;
728 <    e1next = (e1+1)%3;
729 <    e1prev = (e1+2)%3;
730 <    verts[e] = T_NTH_V(SM_NTH_TRI(sm,t_id),e);
731 <    verts[enext] = T_NTH_V(SM_NTH_TRI(sm,t_id),enext);
732 <    verts[eprev] = T_NTH_V(SM_NTH_TRI(sm,t1_id),e1);
733 <    ta_id = smAdd_tri(sm,verts[0],verts[1],verts[2]);
734 <    *add_ptr = push_data(*add_ptr,ta_id);
735 <    verts[e1] = T_NTH_V(SM_NTH_TRI(sm,t1_id),e1);
736 <    verts[e1next] = T_NTH_V(SM_NTH_TRI(sm,t1_id),e1next);
737 <    verts[e1prev] = T_NTH_V(SM_NTH_TRI(sm,t_id),e);
738 <    tb_id = smAdd_tri(sm,verts[0],verts[1],verts[2]);
729 >  int e1,topp_id;
730 >  S_ID p_id;
731 >  TRI *topp;
732 >  FVECT p;
733  
734 <    *add_ptr = push_data(*add_ptr,tb_id);
741 <    ta = SM_NTH_TRI(sm,ta_id);
742 <    tb = SM_NTH_TRI(sm,tb_id);
743 <    t = SM_NTH_TRI(sm,t_id);
744 <    t1 = SM_NTH_TRI(sm,t1_id);
745 <    /* set the neighbors */
746 <    T_NTH_NBR(ta,e) = T_NTH_NBR(t1,e1next);
747 <    T_NTH_NBR(tb,e1) = T_NTH_NBR(t,enext);
748 <    T_NTH_NBR(ta,enext)= tb_id;
749 <    T_NTH_NBR(tb,e1next)= ta_id;
750 <    T_NTH_NBR(ta,eprev)=T_NTH_NBR(t,eprev);
751 <    T_NTH_NBR(tb,e1prev)=T_NTH_NBR(t1,e1prev);    
752 <
753 <  
734 >  topp_id = T_NTH_NBR(t,0);
735   #ifdef DEBUG
736 <        if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,0))) ||
737 <                       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,1))) ||
738 <                       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,2))) ||
739 <                       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,0))) ||
740 <                       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,1))) ||
760 <                       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,2))))
761 <          goto Ltri_error;
762 <        /*
763 <        if(SM_NTH_TRI(sm,T_NTH_NBR(ta,0))==SM_NTH_TRI(sm,T_NTH_NBR(ta,1)) ||
764 <           SM_NTH_TRI(sm,T_NTH_NBR(ta,0))==SM_NTH_TRI(sm,T_NTH_NBR(ta,2)) ||
765 <           SM_NTH_TRI(sm,T_NTH_NBR(ta,1))==SM_NTH_TRI(sm,T_NTH_NBR(ta,2)) ||
766 <           SM_NTH_TRI(sm,T_NTH_NBR(tb,0))==SM_NTH_TRI(sm,T_NTH_NBR(tb,1)) ||
767 <           SM_NTH_TRI(sm,T_NTH_NBR(tb,0))==SM_NTH_TRI(sm,T_NTH_NBR(tb,2)) ||
768 <           SM_NTH_TRI(sm,T_NTH_NBR(tb,1))==SM_NTH_TRI(sm,T_NTH_NBR(tb,2)) )
769 <          goto Ltri_error;
770 <          */
736 >  if(topp_id==INVALID)
737 >  {
738 >    eputs("Invalid tri id smRestore_delaunay()\n");
739 >    return;
740 >  }
741   #endif
742 <    /* Reset neighbor pointers of original neighbors */
773 <    n = SM_NTH_TRI(sm,T_NTH_NBR(t,enext));
742 >  topp = SM_NTH_TRI(sm,topp_id);
743   #ifdef DEBUG
744 <        if(!T_IS_VALID(n))
745 <          goto Ltri_error;
744 >  if(!T_IS_VALID(topp))
745 >  {
746 >    eputs("Invalid tri smRestore_delaunay()\n");
747 >    return;
748 >  }
749   #endif
750 +  e1 = T_NTH_NBR_PTR(t_id,topp);
751 +  p_id = T_NTH_V(topp,e1);
752 +
753 +  VSUB(p,SM_NTH_WV(sm,p_id),SM_VIEW_CENTER(sm));
754 +  normalize(p);
755 +  if(pt_in_cone(p,c,b,a))
756 +  {
757 +    int e1next,e1prev;
758 +    int ta_id,tb_id;
759 +    TRI *ta,*tb,*n;
760 +    e1next = (e1 + 1)%3;
761 +    e1prev = (e1 + 2)%3;
762 +    ta_id = smAdd_tri(sm,a_id,b_id,p_id,&ta);
763 +    tb_id = smAdd_tri(sm,a_id,p_id,c_id,&tb);
764 +
765 +    T_NTH_NBR(ta,0) = T_NTH_NBR(topp,e1next);
766 +    T_NTH_NBR(ta,1) = tb_id;
767 +    T_NTH_NBR(ta,2) = T_NTH_NBR(t,2);
768 +  
769 +    T_NTH_NBR(tb,0) = T_NTH_NBR(topp,e1prev);
770 +    T_NTH_NBR(tb,1) = T_NTH_NBR(t,1);
771 +    T_NTH_NBR(tb,2) = ta_id;
772 +  
773 +    n = SM_NTH_TRI(sm,T_NTH_NBR(t,1));
774      T_NTH_NBR(n,T_NTH_NBR_PTR(t_id,n)) = tb_id;
775 <    n = SM_NTH_TRI(sm,T_NTH_NBR(t,eprev));
780 < #ifdef DEBUG
781 <        if(!T_IS_VALID(n))
782 <          goto Ltri_error;
783 < #endif
775 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t,2));
776      T_NTH_NBR(n,T_NTH_NBR_PTR(t_id,n)) = ta_id;
777 +    n = SM_NTH_TRI(sm,T_NTH_NBR(topp,e1next));
778 +    T_NTH_NBR(n,T_NTH_NBR_PTR(topp_id,n)) = ta_id;
779 +    n = SM_NTH_TRI(sm,T_NTH_NBR(topp,e1prev));
780 +    T_NTH_NBR(n,T_NTH_NBR_PTR(topp_id,n)) = tb_id;
781 +    
782 +    smDelete_tri(sm,t_id,t);
783 +    smDelete_tri(sm,topp_id,topp);
784  
786    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,e1next));
785   #ifdef DEBUG
786 <        if(!T_IS_VALID(n))
787 <          goto Ltri_error;
788 < #endif
789 <    T_NTH_NBR(n,T_NTH_NBR_PTR(t1_id,n)) = ta_id;
790 <    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,e1prev));
791 < #ifdef DEBUG
792 <        if(!T_IS_VALID(n))
793 <          goto Ltri_error;
794 < #endif
795 <    T_NTH_NBR(n,T_NTH_NBR_PTR(t1_id,n)) = tb_id;
786 > #if DEBUG > 1
787 >    if(T_NTH_NBR(ta,0) == T_NTH_NBR(ta,1)  ||
788 >       T_NTH_NBR(ta,1) == T_NTH_NBR(ta,2)  ||
789 >       T_NTH_NBR(ta,0) == T_NTH_NBR(ta,2))
790 >      error(CONSISTENCY,"Invalid topology\n");
791 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,0))) ||
792 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,1))) ||
793 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(ta,2))))
794 >      error(CONSISTENCY,"Invalid topology\n");
795 >    n = SM_NTH_TRI(sm,T_NTH_NBR(ta,0));
796 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
797 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
798 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
799 >      error(CONSISTENCY,"Invalid topology\n");
800 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
801 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
802 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
803 >      error(CONSISTENCY,"Invalid topology\n");
804 >    n = SM_NTH_TRI(sm,T_NTH_NBR(ta,1));
805 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
806 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
807 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
808 >      error(CONSISTENCY,"Invalid topology\n");
809 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
810 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
811 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
812 >      error(CONSISTENCY,"Invalid topology\n");
813 >    n = SM_NTH_TRI(sm,T_NTH_NBR(ta,2));
814 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
815 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
816 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
817 >      error(CONSISTENCY,"Invalid topology\n");
818 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
819 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
820 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
821 >      error(CONSISTENCY,"Invalid topology\n");
822 >    n = SM_NTH_TRI(sm,T_NTH_NBR(tb,0));
823 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
824 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
825 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
826 >      error(CONSISTENCY,"Invalid topology\n");
827 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
828 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
829 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
830 >      error(CONSISTENCY,"Invalid topology\n");
831 >    n = SM_NTH_TRI(sm,T_NTH_NBR(tb,1));
832 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
833 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
834 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
835 >      error(CONSISTENCY,"Invalid topology\n");
836 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
837 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
838 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
839 >      error(CONSISTENCY,"Invalid topology\n");
840 >    n = SM_NTH_TRI(sm,T_NTH_NBR(tb,2));
841 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
842 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
843 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
844 >      error(CONSISTENCY,"Invalid topology\n");
845 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
846 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
847 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
848 >      error(CONSISTENCY,"Invalid topology\n");
849  
850 < #ifdef DEBUG
851 <    T_VALID_FLAG(SM_NTH_TRI(sm,t_id))=-1;
852 <    T_VALID_FLAG(SM_NTH_TRI(sm,t1_id))=-1;
853 < #endif
850 >    if(T_NTH_NBR(tb,0) == T_NTH_NBR(tb,1)  ||
851 >       T_NTH_NBR(tb,1) == T_NTH_NBR(tb,2)  ||
852 >       T_NTH_NBR(tb,0) == T_NTH_NBR(tb,2))
853 >      error(CONSISTENCY,"Invalid topology\n");
854 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,0))) ||
855 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,1))) ||
856 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(tb,2))))
857 >      error(CONSISTENCY,"Invalid topology\n");
858  
859 <    remove_from_list(t_id,add_ptr);
860 <    remove_from_list(t1_id,add_ptr);
861 <
862 <    smDelete_tri(sm,t_id);
808 <    smDelete_tri(sm,t1_id);
809 <
810 <    *tn_id = ta_id;
811 <    *tn1_id = tb_id;
812 <    
813 < #ifdef DEBUG
814 <    if(T_NTH_NBR(SM_NTH_TRI(sm,ta_id),0)== -1 ||
815 <       T_NTH_NBR(SM_NTH_TRI(sm,ta_id),1)== -1 ||
816 <       T_NTH_NBR(SM_NTH_TRI(sm,ta_id),2)== -1 ||
817 <       T_NTH_NBR(SM_NTH_TRI(sm,tb_id),0)== -1 ||
818 <       T_NTH_NBR(SM_NTH_TRI(sm,tb_id),1)== -1 ||
819 <       T_NTH_NBR(SM_NTH_TRI(sm,tb_id),2)== -1)
820 <      goto Ltri_error;
821 <  if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,ta_id),0))) ||
822 <     !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,ta_id),1))) ||
823 <     !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,ta_id),2))))
824 <      goto Ltri_error;
825 <  if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,tb_id),0))) ||
826 <     !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,tb_id),1))) ||
827 <     !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,tb_id),2))))
828 <      goto Ltri_error;
829 < #endif
830 <    return;
831 < Ltri_error:
832 <    error(CONSISTENCY,"Invalid tri: smTris_swap_edge()\n");
833 < }
834 <
835 < smUpdate_locator(sm,add_list)
836 < SM *sm;
837 < LIST *add_list;
838 < {
839 <  int t_id;
840 <  TRI *t;
841 <  
842 <  while(add_list)
843 <  {
844 <    t_id = pop_list(&add_list);
845 <    t = SM_NTH_TRI(sm,t_id);
846 <    smLocator_add_tri(sm,t_id,T_NTH_V(t,0),T_NTH_V(t,1),T_NTH_V(t,2));
859 > #endif  
860 > #endif  
861 >    smRestore_Delaunay(sm,a,b,p,ta,ta_id,a_id,b_id,p_id);
862 >    smRestore_Delaunay(sm,a,p,c,tb,tb_id,a_id,p_id,c_id);
863    }
864   }
865  
850 int
851 smFix_tris(sm,id,tlist,add_list)
852 SM *sm;
853 int id;
854 LIST *tlist,*add_list;
855 {
856    TRI *t,*t_opp;
857    FVECT p,p0,p1,p2;
858    int e,e1,swapped = 0;
859    int t_id,t_opp_id;
860    LIST *del_list=NULL,*lptr;
866  
862    VSUB(p,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
863    while(tlist)
864    {
865        t_id = pop_list(&tlist);
866 #ifdef DEBUG
867        if(t_id==INVALID || t_id > smMesh->num_tri)
868          error(CONSISTENCY,"Invalid tri id smFix_tris()\n");
869 #endif
870        t = SM_NTH_TRI(sm,t_id);
871        e = T_WHICH_V(t,id);
872        t_opp_id = T_NTH_NBR(t,e);
873 #ifdef DEBUG
874        if(t_opp_id==INVALID || t_opp_id > smMesh->num_tri)
875          error(CONSISTENCY,"Invalid tri id smFix_tris()\n");
876 #endif
877        t_opp = SM_NTH_TRI(sm,t_opp_id);
878 #ifdef DEBUG
879        if(!T_IS_VALID(t_opp))
880          error(CONSISTENCY,"Invalid trismFix_tris()\n");
881 #endif
882        smDir(sm,p0,T_NTH_V(t_opp,0));
883        smDir(sm,p1,T_NTH_V(t_opp,1));
884        smDir(sm,p2,T_NTH_V(t_opp,2));
885        if(point_in_cone(p,p0,p1,p2))
886        {
887            swapped = 1;
888            e1 = T_NTH_NBR_PTR(t_id,t_opp);
889            /* check list for t_opp and Remove if there */
890            remove_from_list(t_opp_id,&tlist);
891            smTris_swap_edge(sm,t_id,t_opp_id,e,e1,&t_id,&t_opp_id,
892                             &add_list);
893            tlist = push_data(tlist,t_id);
894            tlist = push_data(tlist,t_opp_id);
895        }
896    }
897
898    smUpdate_locator(sm,add_list);
899
900    return(swapped);
901 }
902
867   /* Give the vertex "id" and a triangle "t" that it belongs to- return the
868     next nbr in a counter clockwise order about vertex "id"
869   */
# Line 907 | Line 871 | int
871   smTri_next_ccw_nbr(sm,t,id)
872   SM *sm;
873   TRI *t;
874 < int id;
874 > S_ID id;
875   {
876    int which;
877    int nbr_id;
# Line 934 | Line 898 | int id;
898  
899   /* Locate the point-id in the point location structure: */
900   int
901 < smInsert_samp(sm,s_id,tri_id,on,which)
901 > smInsert_samp_mesh(sm,s_id,tri_id,a,b,c,d,on,which)
902     SM *sm;
903 <   int s_id,tri_id;
903 >   S_ID s_id;
904 >   int tri_id;
905 >   FVECT a,b,c,d;
906     int on,which;
907   {
908 <    int v_id[3],topp_id,i;
909 <    int t0_id,t1_id,t2_id,t3_id,replaced;
910 <    int e0,e1,e2,opp0,opp1,opp2,opp_id;
911 <    LIST *tlist,*add_list;
912 <    FVECT npt;
913 <    TRI *tri,*nbr,*topp;
908 >    S_ID v_id[3],opp_id;
909 >    int topp_id,i;
910 >    int t0_id,t1_id,t2_id,t3_id;
911 >    int e0,e1,e2,opp0,opp1,opp2;
912 >    TRI *tri,*nbr,*topp,*t0,*t1,*t2,*t3;
913 >    FVECT e;
914  
949    add_list = NULL;
915      for(i=0; i<3;i++)
916        v_id[i] = T_NTH_V(SM_NTH_TRI(sm,tri_id),i);
917  
918      /* If falls on existing edge */
919      if(on == ON_E)
920      {
921 +      FVECT n,vopp;
922 +      double dp;
923 +
924        tri = SM_NTH_TRI(sm,tri_id);
925        topp_id = T_NTH_NBR(tri,which);
926        e0 = which;
# Line 960 | Line 928 | smInsert_samp(sm,s_id,tri_id,on,which)
928        e2 = (which+2)%3;
929        topp = SM_NTH_TRI(sm,topp_id);
930        opp0 = T_NTH_NBR_PTR(tri_id,topp);
963      opp_id = T_NTH_V(topp,opp0);
931        opp1 = (opp0+1)%3;
932        opp2 = (opp0+2)%3;
933 +
934 +      opp_id = T_NTH_V(topp,opp0);
935        
936        /* Create 4 triangles */
937        /*      /|\             /|\
# Line 972 | Line 941 | smInsert_samp(sm,s_id,tri_id,on,which)
941               \ | /           \ | /
942                \|/             \|/
943       */
944 <        t0_id = smAdd_tri(sm,s_id,v_id[e0],v_id[e1]);
945 <        add_list = push_data(add_list,t0_id);
946 <        t1_id = smAdd_tri(sm,s_id,v_id[e2],v_id[e0]);
947 <        add_list = push_data(add_list,t1_id);
979 <        t2_id = smAdd_tri(sm,s_id,v_id[e1],opp_id);
980 <        add_list = push_data(add_list,t2_id);
981 <        t3_id = smAdd_tri(sm,s_id,opp_id,v_id[e2]);
982 <        add_list = push_data(add_list,t3_id);
983 <        /* CAUTION: tri must be assigned after Add_tri: pointers may change*/
984 <        tri = SM_NTH_TRI(sm,tri_id);
985 <        topp =SM_NTH_TRI(sm,topp_id);
944 >      t0_id = smAdd_tri(sm,s_id,v_id[e0],v_id[e1],&t0);
945 >      t1_id = smAdd_tri(sm,s_id,v_id[e2],v_id[e0],&t1);
946 >      t2_id = smAdd_tri(sm,s_id,v_id[e1],opp_id,&t2);
947 >      t3_id = smAdd_tri(sm,s_id,opp_id,v_id[e2],&t3);
948  
949 <        /* Set the neighbor pointers for the new tris */
950 <        T_NTH_NBR(SM_NTH_TRI(sm,t0_id),0) = T_NTH_NBR(tri,e2);
951 <        T_NTH_NBR(SM_NTH_TRI(sm,t0_id),1) = t2_id;
952 <        T_NTH_NBR(SM_NTH_TRI(sm,t0_id),2) = t1_id;
953 <        T_NTH_NBR(SM_NTH_TRI(sm,t1_id),0) = T_NTH_NBR(tri,e1);
954 <        T_NTH_NBR(SM_NTH_TRI(sm,t1_id),1) = t0_id;
955 <        T_NTH_NBR(SM_NTH_TRI(sm,t1_id),2) = t3_id;
956 <        T_NTH_NBR(SM_NTH_TRI(sm,t2_id),0) = T_NTH_NBR(topp,opp1);
957 <        T_NTH_NBR(SM_NTH_TRI(sm,t2_id),1) = t3_id;
958 <        T_NTH_NBR(SM_NTH_TRI(sm,t2_id),2) = t0_id;
959 <        T_NTH_NBR(SM_NTH_TRI(sm,t3_id),0) = T_NTH_NBR(topp,opp2);
960 <        T_NTH_NBR(SM_NTH_TRI(sm,t3_id),1) = t1_id;
961 <        T_NTH_NBR(SM_NTH_TRI(sm,t3_id),2) = t2_id;
949 >      /* Set the neighbor pointers for the new tris */
950 >      T_NTH_NBR(t0,0) = T_NTH_NBR(tri,e2);
951 >      T_NTH_NBR(t0,1) = t2_id;
952 >      T_NTH_NBR(t0,2) = t1_id;
953 >      T_NTH_NBR(t1,0) = T_NTH_NBR(tri,e1);
954 >      T_NTH_NBR(t1,1) = t0_id;
955 >      T_NTH_NBR(t1,2) = t3_id;
956 >      T_NTH_NBR(t2,0) = T_NTH_NBR(topp,opp1);
957 >      T_NTH_NBR(t2,1) = t3_id;
958 >      T_NTH_NBR(t2,2) = t0_id;
959 >      T_NTH_NBR(t3,0) = T_NTH_NBR(topp,opp2);
960 >      T_NTH_NBR(t3,1) = t1_id;
961 >      T_NTH_NBR(t3,2) = t2_id;
962  
963          /* Reset the neigbor pointers for the neighbors of the original */
964          nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,e1));
# Line 1007 | Line 969 | smInsert_samp(sm,s_id,tri_id,on,which)
969          T_NTH_NBR(nbr,T_NTH_NBR_PTR(topp_id,nbr)) = t2_id;
970          nbr = SM_NTH_TRI(sm,T_NTH_NBR(topp,opp2));
971          T_NTH_NBR(nbr,T_NTH_NBR_PTR(topp_id,nbr)) = t3_id;
972 +
973   #ifdef DEBUG
974 <    T_VALID_FLAG(SM_NTH_TRI(sm,tri_id))=-1;
975 <    T_VALID_FLAG(SM_NTH_TRI(sm,topp_id))=-1;
976 < #endif
977 <        tlist = push_data(NULL,t0_id);
978 <        tlist = push_data(tlist,t1_id);
979 <        tlist = push_data(tlist,t2_id);
980 <        tlist = push_data(tlist,t3_id);
974 > #if DEBUG > 1
975 > {
976 >   TRI *n;
977 >
978 >    if(T_NTH_NBR(t0,0) == T_NTH_NBR(t0,1)  ||
979 >       T_NTH_NBR(t0,1) == T_NTH_NBR(t0,2)  ||
980 >       T_NTH_NBR(t0,0) == T_NTH_NBR(t0,2))
981 >      error(CONSISTENCY,"Invalid topology\n");
982 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,0))) ||
983 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,1))) ||
984 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,2))))
985 >      error(CONSISTENCY,"Invalid topology\n");
986 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,0));
987 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
988 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
989 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
990 >      error(CONSISTENCY,"Invalid topology\n");
991 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
992 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
993 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
994 >      error(CONSISTENCY,"Invalid topology\n");
995 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,1));
996 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
997 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
998 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
999 >      error(CONSISTENCY,"Invalid topology\n");
1000 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1001 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1002 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1003 >      error(CONSISTENCY,"Invalid topology\n");
1004 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,2));
1005 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1006 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1007 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1008 >      error(CONSISTENCY,"Invalid topology\n");
1009 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1010 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1011 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1012 >      error(CONSISTENCY,"Invalid topology\n");
1013 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,0));
1014 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1015 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1016 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1017 >      error(CONSISTENCY,"Invalid topology\n");
1018 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1019 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1020 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1021 >      error(CONSISTENCY,"Invalid topology\n");
1022 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,1));
1023 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1024 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1025 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1026 >      error(CONSISTENCY,"Invalid topology\n");
1027 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1028 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1029 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1030 >      error(CONSISTENCY,"Invalid topology\n");
1031 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,2));
1032 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1033 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1034 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1035 >      error(CONSISTENCY,"Invalid topology\n");
1036 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1037 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1038 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1039 >      error(CONSISTENCY,"Invalid topology\n");
1040 >
1041 >    if(T_NTH_NBR(t1,0) == T_NTH_NBR(t1,1)  ||
1042 >       T_NTH_NBR(t1,1) == T_NTH_NBR(t1,2)  ||
1043 >       T_NTH_NBR(t1,0) == T_NTH_NBR(t1,2))
1044 >      error(CONSISTENCY,"Invalid topology\n");
1045 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,0))) ||
1046 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,1))) ||
1047 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,2))))
1048 >      error(CONSISTENCY,"Invalid topology\n");
1049 >
1050 >    if(T_NTH_NBR(t2,0) == T_NTH_NBR(t2,1)  ||
1051 >       T_NTH_NBR(t2,1) == T_NTH_NBR(t2,2)  ||
1052 >       T_NTH_NBR(t2,0) == T_NTH_NBR(t2,2))
1053 >      error(CONSISTENCY,"Invalid topology\n");
1054 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t2,0))) ||
1055 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t2,1))) ||
1056 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t2,2))))
1057 >      error(CONSISTENCY,"Invalid topology\n");
1058 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t2,0));
1059 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1060 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1061 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1062 >      error(CONSISTENCY,"Invalid topology\n");
1063 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1064 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1065 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1066 >      error(CONSISTENCY,"Invalid topology\n");
1067 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t2,1));
1068 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1069 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1070 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1071 >      error(CONSISTENCY,"Invalid topology\n");
1072 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1073 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1074 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1075 >      error(CONSISTENCY,"Invalid topology\n");
1076 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t2,2));
1077 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1078 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1079 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1080 >      error(CONSISTENCY,"Invalid topology\n");
1081 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1082 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1083 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1084 >      error(CONSISTENCY,"Invalid topology\n");
1085 >
1086 >    if(T_NTH_NBR(t3,0) == T_NTH_NBR(t3,1)  ||
1087 >       T_NTH_NBR(t3,1) == T_NTH_NBR(t3,2)  ||
1088 >       T_NTH_NBR(t3,0) == T_NTH_NBR(t3,2))
1089 >      error(CONSISTENCY,"Invalid topology\n");
1090 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t3,0))) ||
1091 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t3,1))) ||
1092 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t3,2))))
1093 >      error(CONSISTENCY,"Invalid topology\n");
1094 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t3,0));
1095 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1096 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1097 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1098 >      error(CONSISTENCY,"Invalid topology\n");
1099 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1100 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1101 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1102 >      error(CONSISTENCY,"Invalid topology\n");
1103 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t3,1));
1104 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1105 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1106 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1107 >      error(CONSISTENCY,"Invalid topology\n");
1108 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1109 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1110 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1111 >      error(CONSISTENCY,"Invalid topology\n");
1112 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t3,2));
1113 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1114 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1115 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1116 >      error(CONSISTENCY,"Invalid topology\n");
1117 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1118 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1119 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1120 >      error(CONSISTENCY,"Invalid topology\n");
1121 > }
1122 > #endif  
1123 > #endif  
1124 >        smDir(sm,e,opp_id);
1125 >        if(which == 0)
1126 >        {
1127 >          smRestore_Delaunay(sm,a,b,c,t0,t0_id,s_id,v_id[e0],v_id[e1]);
1128 >          smRestore_Delaunay(sm,a,d,b,t1,t1_id,s_id,v_id[e2],v_id[e0]);
1129 >          smRestore_Delaunay(sm,a,c,e,t2,t2_id,s_id,v_id[e1],opp_id);
1130 >          smRestore_Delaunay(sm,a,e,d,t3,t3_id,s_id,opp_id,v_id[e2]);
1131 >        }
1132 >        else
1133 >          if( which==1 )
1134 >          {
1135 >            smRestore_Delaunay(sm,a,c,d,t0,t0_id,s_id,v_id[e0],v_id[e1]);
1136 >            smRestore_Delaunay(sm,a,b,c,t1,t1_id,s_id,v_id[e2],v_id[e0]);
1137 >            smRestore_Delaunay(sm,a,d,e,t2,t2_id,s_id,v_id[e1],opp_id);
1138 >            smRestore_Delaunay(sm,a,e,b,t3,t3_id,s_id,opp_id,v_id[e2]);
1139 >          }
1140 >          else
1141 >          {
1142 >            smRestore_Delaunay(sm,a,d,b,t0,t0_id,s_id,v_id[e0],v_id[e1]);
1143 >            smRestore_Delaunay(sm,a,c,d,t1,t1_id,s_id,v_id[e2],v_id[e0]);
1144 >            smRestore_Delaunay(sm,a,b,e,t2,t2_id,s_id,v_id[e1],opp_id);
1145 >            smRestore_Delaunay(sm,a,e,c,t3,t3_id,s_id,opp_id,v_id[e2]);
1146 >          }
1147 >        smDelete_tri(sm,tri_id,tri);
1148 >        smDelete_tri(sm,topp_id,topp);
1149 >        return(TRUE);
1150      }
1151 <    else
1152 <      {
1153 <        /* Create 3 triangles */
1154 <      /*      / \             /|\
1155 <             /   \           / | \
1156 <            /  *  \  ---->  /  |  \
1157 <           /       \       /  / \  \
1158 <          /         \     / /     \ \
1159 <         /___________\   //_________\\
1160 <      */
1151 > Linsert_in_tri:
1152 >    /* Create 3 triangles */
1153 >    /*      / \             /|\
1154 >           /   \           / | \
1155 >          /  *  \  ---->  /  |  \
1156 >         /       \       /  / \  \
1157 >        /         \       / /     \ \
1158 >        ___________\   //_________\\
1159 >    */
1160 >    tri = SM_NTH_TRI(sm,tri_id);
1161  
1162 <        t0_id = smAdd_tri(sm,s_id,v_id[0],v_id[1]);
1163 <        add_list = push_data(add_list,t0_id);
1164 <        t1_id = smAdd_tri(sm,s_id,v_id[1],v_id[2]);
1033 <        add_list = push_data(add_list,t1_id);
1034 <        t2_id = smAdd_tri(sm,s_id,v_id[2],v_id[0]);
1035 <        add_list = push_data(add_list,t2_id);
1162 >    t0_id = smAdd_tri(sm,s_id,v_id[0],v_id[1],&t0);
1163 >    t1_id = smAdd_tri(sm,s_id,v_id[1],v_id[2],&t1);
1164 >    t2_id = smAdd_tri(sm,s_id,v_id[2],v_id[0],&t2);
1165  
1166 <        /* CAUTION: tri must be assigned after Add_tri: pointers may change*/
1167 <        tri = SM_NTH_TRI(sm,tri_id);
1168 <        /* Set the neighbor pointers for the new tris */
1169 <        T_NTH_NBR(SM_NTH_TRI(sm,t0_id),0) = T_NTH_NBR(tri,2);
1170 <        T_NTH_NBR(SM_NTH_TRI(sm,t0_id),1) = t1_id;
1171 <        T_NTH_NBR(SM_NTH_TRI(sm,t0_id),2) = t2_id;
1172 <        T_NTH_NBR(SM_NTH_TRI(sm,t1_id),0) = T_NTH_NBR(tri,0);
1173 <        T_NTH_NBR(SM_NTH_TRI(sm,t1_id),1) = t2_id;
1174 <        T_NTH_NBR(SM_NTH_TRI(sm,t1_id),2) = t0_id;
1175 <        T_NTH_NBR(SM_NTH_TRI(sm,t2_id),0) = T_NTH_NBR(tri,1);
1047 <        T_NTH_NBR(SM_NTH_TRI(sm,t2_id),1) = t0_id;
1048 <        T_NTH_NBR(SM_NTH_TRI(sm,t2_id),2) = t1_id;
1166 >    /* Set the neighbor pointers for the new tris */
1167 >    T_NTH_NBR(t0,0) = T_NTH_NBR(tri,2);
1168 >    T_NTH_NBR(t0,1) = t1_id;
1169 >    T_NTH_NBR(t0,2) = t2_id;
1170 >    T_NTH_NBR(t1,0) = T_NTH_NBR(tri,0);
1171 >    T_NTH_NBR(t1,1) = t2_id;
1172 >    T_NTH_NBR(t1,2) = t0_id;
1173 >    T_NTH_NBR(t2,0) = T_NTH_NBR(tri,1);
1174 >    T_NTH_NBR(t2,1) = t0_id;
1175 >    T_NTH_NBR(t2,2) = t1_id;
1176          
1177 <        /* Reset the neigbor pointers for the neighbors of the original */
1178 <        nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,0));
1179 <        T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t1_id;
1180 <        nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,1));
1181 <        T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t2_id;
1182 <        nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,2));
1183 <        T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t0_id;
1177 >    /* Reset the neigbor pointers for the neighbors of the original */
1178 >    nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,0));
1179 >    T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t1_id;
1180 >    nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,1));
1181 >    T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t2_id;
1182 >    nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,2));
1183 >    T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t0_id;
1184 > #ifdef DEBUG
1185 > #if DEBUG > 1
1186 > {
1187 >  TRI *n;
1188 >    if(T_NTH_NBR(t0,0) == T_NTH_NBR(t0,1)  ||
1189 >       T_NTH_NBR(t0,1) == T_NTH_NBR(t0,2)  ||
1190 >       T_NTH_NBR(t0,0) == T_NTH_NBR(t0,2))
1191 >      error(CONSISTENCY,"Invalid topology\n");
1192 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,0))) ||
1193 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,1))) ||
1194 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t0,2))))
1195 >      error(CONSISTENCY,"Invalid topology\n");
1196 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,0));
1197 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1198 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1199 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1200 >      error(CONSISTENCY,"Invalid topology\n");
1201 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1202 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1203 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1204 >      error(CONSISTENCY,"Invalid topology\n");
1205 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,1));
1206 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1207 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1208 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1209 >      error(CONSISTENCY,"Invalid topology\n");
1210 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1211 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1212 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1213 >      error(CONSISTENCY,"Invalid topology\n");
1214 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t0,2));
1215 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1216 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1217 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1218 >      error(CONSISTENCY,"Invalid topology\n");
1219 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1220 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1221 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1222 >      error(CONSISTENCY,"Invalid topology\n");
1223 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,0));
1224 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1225 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1226 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1227 >      error(CONSISTENCY,"Invalid topology\n");
1228 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1229 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1230 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1231 >      error(CONSISTENCY,"Invalid topology\n");
1232 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,1));
1233 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1234 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1235 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1236 >      error(CONSISTENCY,"Invalid topology\n");
1237 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1238 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1239 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1240 >      error(CONSISTENCY,"Invalid topology\n");
1241 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t1,2));
1242 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1243 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1244 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1245 >      error(CONSISTENCY,"Invalid topology\n");
1246 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1247 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1248 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1249 >      error(CONSISTENCY,"Invalid topology\n");
1250  
1251 <        tlist = push_data(NULL,t0_id);
1252 <        tlist = push_data(tlist,t1_id);
1253 <        tlist = push_data(tlist,t2_id);
1254 <      }
1255 <    /* Fix up the new triangles*/
1256 <    smFix_tris(sm,s_id,tlist,add_list);
1251 >    if(T_NTH_NBR(t1,0) == T_NTH_NBR(t1,1)  ||
1252 >       T_NTH_NBR(t1,1) == T_NTH_NBR(t1,2)  ||
1253 >       T_NTH_NBR(t1,0) == T_NTH_NBR(t1,2))
1254 >      error(CONSISTENCY,"Invalid topology\n");
1255 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,0))) ||
1256 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,1))) ||
1257 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t1,2))))
1258 >      error(CONSISTENCY,"Invalid topology\n");
1259  
1260 <    smDelete_tri(sm,tri_id);
1261 <    if(on==ON_E)
1262 <      smDelete_tri(sm,topp_id);
1260 >    if(T_NTH_NBR(t2,0) == T_NTH_NBR(t2,1)  ||
1261 >       T_NTH_NBR(t2,1) == T_NTH_NBR(t2,2)  ||
1262 >       T_NTH_NBR(t2,0) == T_NTH_NBR(t2,2))
1263 >      error(CONSISTENCY,"Invalid topology\n");
1264 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t2,0))) ||
1265 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t2,1))) ||
1266 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(t2,2))))
1267 >      error(CONSISTENCY,"Invalid topology\n");
1268 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t2,0));
1269 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1270 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1271 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1272 >      error(CONSISTENCY,"Invalid topology\n");
1273 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1274 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1275 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1276 >      error(CONSISTENCY,"Invalid topology\n");
1277 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t2,1));
1278 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1279 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1280 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1281 >      error(CONSISTENCY,"Invalid topology\n");
1282 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1283 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1284 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1285 >      error(CONSISTENCY,"Invalid topology\n");
1286 >    n = SM_NTH_TRI(sm,T_NTH_NBR(t2,2));
1287 >    if(T_NTH_NBR(n,0) == T_NTH_NBR(n,1)  ||
1288 >       T_NTH_NBR(n,1) == T_NTH_NBR(n,2)  ||
1289 >       T_NTH_NBR(n,0) == T_NTH_NBR(n,2))
1290 >      error(CONSISTENCY,"Invalid topology\n");
1291 >    if(!T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,0))) ||
1292 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,1))) ||
1293 >       !T_IS_VALID(SM_NTH_TRI(sm,T_NTH_NBR(n,2))))
1294 >      error(CONSISTENCY,"Invalid topology\n");
1295 > }
1296 > #endif
1297 > #endif
1298 >  smRestore_Delaunay(sm,a,b,c,t0,t0_id,s_id,v_id[0],v_id[1]);
1299 >  smRestore_Delaunay(sm,a,c,d,t1,t1_id,s_id,v_id[1],v_id[2]);
1300 >  smRestore_Delaunay(sm,a,d,b,t2,t2_id,s_id,v_id[2],v_id[0]);
1301 >  smDelete_tri(sm,tri_id,tri);
1302 >  return(TRUE);
1303  
1069    return(TRUE);
1070
1304   #ifdef DEBUG
1305   Ladderror:
1306 <        error(CONSISTENCY,"Invalid tri: smInsert_samp()\n");
1306 >        error(CONSISTENCY,"Invalid tri: smInsert_samp_mesh()\n");
1307   #endif
1308   }    
1309  
1310 +
1311   int
1312 < smTri_in_set(sm,p,optr,onptr,whichptr)
1312 > smWalk(sm,p,t_id,t,onptr,wptr,from,on_edge,a,b)
1313   SM *sm;
1314   FVECT p;
1315 < OBJECT *optr;
1316 < int *onptr,*whichptr;
1315 > int t_id;
1316 > TRI *t;
1317 > int *onptr,*wptr;
1318 > int from;
1319 > double on_edge;
1320 > FVECT a,b;
1321   {
1084  int i,t_id,on0,on1,on2;
1322    FVECT n,v0,v1,v2;
1323 <  TRI *t;
1324 <  double side;
1323 >  TRI *tn;
1324 >  double on0,on1,on2;
1325 >  int tn_id;
1326  
1327 <  for (i = QT_SET_CNT(optr),optr = QT_SET_PTR(optr);i > 0; i--)
1328 <  {
1329 <    /* Find the first triangle that pt falls */
1330 <    t_id = QT_SET_NEXT_ELEM(optr);
1331 <    t = SM_NTH_TRI(sm,t_id);
1332 <    if(!T_IS_VALID(t))
1333 <      continue;
1327 > #ifdef DEBUG
1328 >  Wcnt++;
1329 >  on0 = on1 = on2 = 10;
1330 > #endif
1331 >  switch(from){
1332 >  case 0:
1333 >    on0 = on_edge;
1334 >    /* Havnt calculate v0 yet: check for equality first */
1335      VSUB(v0,SM_T_NTH_WV(sm,t,0),SM_VIEW_CENTER(sm));
1336 <    VSUB(v1,SM_T_NTH_WV(sm,t,1),SM_VIEW_CENTER(sm));
1337 <    VSUB(v2,SM_T_NTH_WV(sm,t,2),SM_VIEW_CENTER(sm));
1336 >    normalize(v0);
1337 >    if(EQUAL_VEC3(v0,p)){ *onptr = ON_V; *wptr = 0; return(t_id); }
1338 >    /* determine what side of edge v0v1 (v0a) point is on*/
1339 >    VCROSS(n,v0,a);
1340 >    normalize(n);
1341 >    on2 = DOT(n,p);
1342 >    if(on2 > 0.0)
1343 >      {
1344 >        tn_id = T_NTH_NBR(t,2);
1345 >        tn = SM_NTH_TRI(sm,tn_id);
1346 >        return(smWalk(sm,p,tn_id,tn,onptr,wptr,T_NTH_NBR_PTR(t_id,tn),
1347 >                      -on2,a,v0));
1348 >      }
1349 >    else
1350 >      if(on2 >= -EDGE_EPS)
1351 >      {
1352 >        if(on0 >= -EDGE_EPS)
1353 >        {  
1354 >          /* Could either be epsilon of vertex 1*/
1355 >        /* dot(p,a) == cos angle between them. Want length of chord x to
1356 >           be <= EDGE_EPS for test. if the perpendicular bisector of pa
1357 >           is d, cos@/2 = d/1, with right triangle d^2 + (x/2)^2 = 1
1358 >           or x^2 = 4(1-cos^2@/2) = 4(1- (1 + cos@)/2) =  2 -2cos@,
1359 >            so test is if  x^2 <= VERT_EPS*VERT_EPS,
1360 >           2 - 2cos@ <= VERT_EPS*VERT_EPS, or... */
1361 >          if(DOT(p,a) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1362 >            {*onptr = ON_P; *wptr = 1;}
1363 >          else
1364 >            if(on2 > on0)
1365 >            {*onptr = ON_E; *wptr = 2;}
1366 >            else
1367 >            {*onptr = ON_E; *wptr = 0;}
1368 >          return(t_id);
1369 >        }
1370 >      }
1371  
1372 <    if(EQUAL_VEC3(v0,p))
1372 >    VCROSS(n,b,v0);
1373 >    normalize(n);
1374 >    on1 = DOT(n,p);
1375 >    if(on1 > 0.0)
1376      {
1377 <      *onptr = ON_V;
1378 <      *whichptr = 0;
1379 <      return(t_id);
1377 >      tn_id = T_NTH_NBR(t,1);
1378 >      tn = SM_NTH_TRI(sm,tn_id);
1379 >      return(smWalk(sm,p,tn_id,tn,onptr,wptr,T_NTH_NBR_PTR(t_id,tn),
1380 >                   -on1,v0,b));
1381      }
1382 <    if(EQUAL_VEC3(v1,p))
1382 >    else
1383 >      if(on1 >= -EDGE_EPS)
1384 >      { /* Either on v0 or on edge 1 */
1385 >        if(on0 >= -EDGE_EPS)
1386 >        {
1387 >        if(DOT(p,b) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1388 >            {*onptr = ON_P; *wptr = 2;}
1389 >        else
1390 >          if(on1 > on0)
1391 >            {*onptr = ON_E; *wptr = 1;}
1392 >          else
1393 >            {*onptr = ON_E; *wptr = 0;}
1394 >          return(t_id);
1395 >        }
1396 >        if(on2 >= -EDGE_EPS)
1397 >        {
1398 >        if(DOT(p,v0) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1399 >            {*onptr = ON_P; *wptr = 0;}
1400 >        else
1401 >          if( on1 > on2)
1402 >            {*onptr = ON_E; *wptr = 1;}
1403 >          else
1404 >            {*onptr = ON_E; *wptr = 2;}
1405 >          return(t_id);
1406 >        }
1407 >        *onptr = ON_E; *wptr = 1;
1408 >        return(t_id);
1409 >      }
1410 >    else
1411 >      if(on2 >= -EDGE_EPS)
1412 >      {
1413 >        *onptr = ON_E; *wptr = 2;
1414 >        return(t_id);
1415 >      }
1416 >      if(on0 >= -EDGE_EPS)
1417 >      {
1418 >        *onptr = ON_E; *wptr = 0;
1419 >        return(t_id);
1420 >      }
1421 >    *onptr = IN_T;
1422 >    return(t_id);
1423 >  case 1:
1424 >    on1 = on_edge;
1425 >    /* Havnt calculate v1 yet: check for equality first */
1426 >    VSUB(v1,SM_T_NTH_WV(sm,t,1),SM_VIEW_CENTER(sm));
1427 >    normalize(v1);
1428 >    if(EQUAL_VEC3(v1,p)){ *onptr = ON_V; *wptr = 1; return(t_id); }
1429 >    /* determine what side of edge v0v1 (v0a) point is on*/
1430 >    VCROSS(n,b,v1);
1431 >    normalize(n);
1432 >    on2 = DOT(n,p);
1433 >    if(on2 > 0.0)
1434      {
1435 <      *onptr = ON_V;
1436 <      *whichptr = 1;
1437 <      return(t_id);
1435 >      tn_id = T_NTH_NBR(t,2);
1436 >      tn = SM_NTH_TRI(sm,tn_id);
1437 >      return(smWalk(sm,p,tn_id,tn,onptr,wptr,T_NTH_NBR_PTR(t_id,tn),
1438 >                   -on2,v1,b));
1439      }
1440 <    if(EQUAL_VEC3(v2,p))
1440 >    if(on2 >= -EDGE_EPS)
1441      {
1442 <      *onptr = ON_V;
1443 <      *whichptr = 2;
1444 <      return(t_id);
1442 >      if(on1 >= -EDGE_EPS)
1443 >      {
1444 >        if(DOT(p,b) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1445 >          { *onptr = ON_P; *wptr = 0;}
1446 >        else
1447 >          if(on2 > on1)
1448 >          { *onptr = ON_E; *wptr = 2;}
1449 >        else
1450 >          { *onptr = ON_E; *wptr = 1;}
1451 >         return(t_id);
1452 >      }
1453      }
1118    
1119    on0 = on1 =on2 = 0;
1120    VCROSS(n,v0,v1);
1121    side = DOT(n,p);
1122    if(ZERO(side))
1123      on2 = TRUE;
1124    else
1125      if(side >0.0)
1126        continue;
1454  
1455 <    VCROSS(n,v1,v2);
1456 <    side = DOT(n,p);
1457 <    if(ZERO(side))
1458 <      on0 = TRUE;
1455 >    VCROSS(n,v1,a);
1456 >    normalize(n);
1457 >    on0 = DOT(n,p);
1458 >    if(on0 >0.0)
1459 >    {
1460 >      tn_id = T_NTH_NBR(t,0);
1461 >      tn = SM_NTH_TRI(sm,tn_id);
1462 >      return(smWalk(sm,p,tn_id,tn,onptr,wptr,T_NTH_NBR_PTR(t_id,tn),
1463 >                 -on0,a,v1));
1464 >    }
1465 >    else
1466 >    if(on0 >= -EDGE_EPS)
1467 >      { /* Either on v0 or on edge 1 */
1468 >       if(on1 >= -EDGE_EPS)
1469 >       {
1470 >        if(DOT(p,a) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1471 >           {*onptr = ON_P; *wptr = 2;}
1472 >         else
1473 >          if(on0 > on1)
1474 >          { *onptr = ON_E; *wptr = 0;}
1475 >        else
1476 >          { *onptr = ON_E; *wptr = 1;}
1477 >         return(t_id);
1478 >       }
1479 >       if(on2 >= -EDGE_EPS)
1480 >       {
1481 >       if(DOT(p,v1) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1482 >           {*onptr = ON_P; *wptr = 1;}
1483 >         else
1484 >          if(on0 > on2)
1485 >          { *onptr = ON_E; *wptr = 0;}
1486 >        else
1487 >          { *onptr = ON_E; *wptr = 2;}
1488 >         return(t_id);
1489 >       }
1490 >       *onptr = ON_E; *wptr = 0;
1491 >       return(t_id);
1492 >     }
1493      else
1494 <      if(side >0.0)
1495 <        continue;
1494 >      if(on2 >= -EDGE_EPS)
1495 >      {
1496 >        *onptr = ON_E; *wptr = 2;
1497 >        return(t_id);
1498 >      }
1499 >      if(on1 >= -EDGE_EPS)
1500 >      {
1501 >        *onptr = ON_E; *wptr = 1;
1502 >        return(t_id);
1503 >      }
1504 >    *onptr = IN_T;
1505 >    return(t_id);
1506 >  case 2:
1507 >    on2 = on_edge;
1508 >    VSUB(v2,SM_T_NTH_WV(sm,t,2),SM_VIEW_CENTER(sm));
1509 >    normalize(v2);
1510 >    if(EQUAL_VEC3(v2,p)){ *onptr = ON_V; *wptr = 2; return(t_id); }
1511  
1512 <    VCROSS(n,v2,v0);
1513 <    side = DOT(n,p);
1514 <    if(ZERO(side))
1515 <      on1 = TRUE;
1516 <    else
1141 <      if(side >0.0)
1142 <        continue;
1143 <    if(on0)
1512 >    /* determine what side of edge v0v1 (v0a) point is on*/
1513 >    VCROSS(n,b,v2);
1514 >    normalize(n);
1515 >    on0 = DOT(n,p);
1516 >    if(on0 > 0.0)
1517      {
1518 <      if(on1)
1518 >      tn_id = T_NTH_NBR(t,0);
1519 >      tn = SM_NTH_TRI(sm,tn_id);
1520 >      return(smWalk(sm,p,tn_id,tn,onptr,wptr,T_NTH_NBR_PTR(t_id,tn),
1521 >                    -on0,v2,b));
1522 >    }
1523 >    else
1524 >    if(on0 >= -EDGE_EPS)
1525        {
1526 <        *onptr = ON_P;
1527 <        *whichptr = 2;
1526 >        if(on2 >= - EDGE_EPS)
1527 >        {
1528 >        if(DOT(p,b) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1529 >           {*onptr = ON_P; *wptr = 1;}
1530 >         else
1531 >           if(on0 > on2)
1532 >             {*onptr = ON_E; *wptr = 0;}
1533 >           else
1534 >             {*onptr = ON_E; *wptr = 2;}
1535 >         return(t_id);
1536 >        }
1537        }
1538 <      else
1539 <        if(on2)
1538 >    VCROSS(n,v2,a);
1539 >    normalize(n);
1540 >    on1 = DOT(n,p);
1541 >    if(on1 >0.0)
1542 >    {
1543 >      tn_id = T_NTH_NBR(t,1);
1544 >      tn = SM_NTH_TRI(sm,tn_id);
1545 >      return(smWalk(sm,p,tn_id,tn,onptr,wptr,T_NTH_NBR_PTR(t_id,tn),
1546 >                   -on1,a,v2));
1547 >    }
1548 >    else
1549 >    if(on1 >= -EDGE_EPS)
1550 >     { /* Either on v0 or on edge 1 */
1551 >       if(on0 >= - EDGE_EPS)
1552         {
1553 <         *onptr = ON_P;
1554 <         *whichptr = 1;
1553 >        if(DOT(p,v2) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1554 >          {*onptr = ON_P;*wptr = 2;}
1555 >        else
1556 >          if(on1 > on0)
1557 >            {*onptr = ON_E;*wptr = 1;}
1558 >          else
1559 >            {*onptr = ON_E;*wptr = 0;}
1560 >         return(t_id);
1561         }
1562 +       if(on2 >= -EDGE_EPS)
1563 +       {
1564 +        if(DOT(p,a) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1565 +          {*onptr = ON_P;*wptr = 0;}
1566          else
1567 <          {
1568 <            *onptr = ON_E;
1569 <            *whichptr = 0;
1570 <          }
1567 >          if(on1 > on2)
1568 >            {*onptr = ON_E;*wptr = 1;}
1569 >          else
1570 >            {*onptr = ON_E;*wptr = 2;}
1571 >         return(t_id);
1572 >       }
1573 >       *onptr = ON_E;*wptr = 1;
1574 >       return(t_id);
1575 >     }
1576 >    else
1577 >      if(on0 >= -EDGE_EPS)
1578 >      {
1579 >        *onptr = ON_E;*wptr = 0;
1580 >        return(t_id);
1581 >      }
1582 >      if(on2 >= -EDGE_EPS)
1583 >      {
1584 >        *onptr = ON_E;*wptr = 2;
1585 >        return(t_id);
1586 >      }
1587 >    *onptr = IN_T;
1588 >    return(t_id);
1589 >  default:
1590 >    /* First time called: havnt tested anything */
1591 >    /* Check against vertices */
1592 >    VSUB(v0,SM_T_NTH_WV(sm,t,0),SM_VIEW_CENTER(sm));
1593 >    normalize(v0);
1594 >    if(EQUAL_VEC3(v0,p)){ *onptr = ON_V; *wptr = 0; return(t_id); }
1595 >    VSUB(v1,SM_T_NTH_WV(sm,t,1),SM_VIEW_CENTER(sm));
1596 >    normalize(v1);
1597 >    if(EQUAL_VEC3(v1,p)){ *onptr = ON_V; *wptr = 1; return(t_id); }
1598 >    VSUB(v2,SM_T_NTH_WV(sm,t,2),SM_VIEW_CENTER(sm));
1599 >    normalize(v2);
1600 >    if(EQUAL_VEC3(v2,p)){ *onptr = ON_V; *wptr = 2; return(t_id); }
1601 >    
1602 >    VCROSS(n,v0,v1);  /* Check against v0v1, or edge 2 */
1603 >    normalize(n);
1604 >    on2 = DOT(n,p);
1605 >    if(on2 > 0.0)
1606 >    {              /* Outside edge: traverse to nbr 2 */
1607 >      tn_id = T_NTH_NBR(t,2);
1608 >      tn = SM_NTH_TRI(sm,tn_id);
1609 >      return(smWalk(sm,p,tn_id,tn,onptr,wptr,T_NTH_NBR_PTR(t_id,tn),
1610 >                 -on2,v1,v0));
1611 >    }
1612 >    else
1613 >    VCROSS(n,v1,v2);  /* Check against v1v2 or edge 0 */
1614 >    normalize(n);
1615 >    on0 = DOT(n,p);
1616 >    if(on0 > 0.0)
1617 >    {                      /* Outside edge: traverse to nbr 0 */
1618 >      tn_id = T_NTH_NBR(t,0);
1619 >      tn = SM_NTH_TRI(sm,tn_id);
1620 >      return(smWalk(sm,p,tn_id,tn,onptr,wptr,T_NTH_NBR_PTR(t_id,tn),
1621 >                    -on0,v2,v1));
1622 >    }
1623 >    else
1624 >      if(on0 >= -EDGE_EPS)
1625 >      {                /* On the line defining edge 0 */
1626 >      if(on2 >= -EDGE_EPS)
1627 >      {              
1628 >        if(DOT(p,v1) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1629 >          {*onptr = ON_P;*wptr = 1;}
1630 >        else
1631 >          if(on0 > on2)
1632 >          {*onptr = ON_E;*wptr = 0;}
1633 >        else
1634 >          {*onptr = ON_E;*wptr = 2;}
1635 >        return(t_id);
1636 >      }
1637 >    }
1638 >    VCROSS(n,v2,v0);  /* Check against v2v0 or edge 1 */
1639 >    normalize(n);
1640 >    on1 = DOT(n,p);
1641 >    if(on1 >0.0)
1642 >    {                            /* Outside edge: traverse to nbr 1 */
1643 >      tn_id = T_NTH_NBR(t,1);
1644 >      tn = SM_NTH_TRI(sm,tn_id);
1645 >      return(smWalk(sm,p,tn_id,tn,onptr,wptr,T_NTH_NBR_PTR(t_id,tn),
1646 >                    -on1,v0,v2));
1647 >    }
1648 >    else
1649 >    if(on1  >= -EDGE_EPS)
1650 >      {               /* On the line defining edge 1 */
1651 >        if(on2 >= -EDGE_EPS)
1652 >        {              
1653 >        if(DOT(p,v0) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1654 >            {*onptr = ON_P;*wptr = 0;}
1655 >          else
1656 >            if(on1 > on2)
1657 >              {*onptr = ON_E;*wptr = 1;}
1658 >            else
1659 >              {*onptr = ON_E;*wptr = 2;}
1660 >          return(t_id);
1661 >        }
1662 >        if(on0 >= -EDGE_EPS)
1663 >        {              
1664 >        if(DOT(p,v2) >= 1.0-VERT_EPS*VERT_EPS/2.0)
1665 >            {*onptr = ON_P;*wptr = 2;}
1666 >          else
1667 >            if(on1 > on0)
1668 >              {*onptr = ON_E;*wptr = 1;}
1669 >            else
1670 >              {*onptr = ON_E;*wptr = 0;}
1671 >          return(t_id);
1672 >       }
1673 >        *onptr = ON_E; /* Only on edge 1 */
1674 >        *wptr = 1;
1675 >        return(t_id);
1676 >      }
1677 >    else
1678 >      if(on2 >= -EDGE_EPS)           /* Is on edge 2 ? */
1679 >      {
1680 >        *onptr = ON_E; *wptr = 2;
1681 >        return(t_id);
1682 >      }
1683 >    if(on0 >= -EDGE_EPS)            /* Is on edge 0 ? */
1684 >    {
1685 >      *onptr = ON_E; *wptr = 0;
1686        return(t_id);
1687      }
1688 <   if(on1)
1688 >    *onptr = IN_T;   /* Must be interior to t */
1689 >    return(t_id);
1690 >  }
1691 > }
1692 >
1693 > #define check_edges(on0,on1,on2,o,w,t_id) \
1694 >   if(on0){if(on1){*o = ON_P;*w = 2;} else if(on2){*o = ON_P;*w = 1;} \
1695 >                   else {*o = ON_E; *w = 0; } return(t_id); } \
1696 >   if(on1){if(on2){*o = ON_P;*w = 0;} else {*o = ON_E;*w = 1;}return(t_id);}\
1697 >   if(on2){*o = ON_E;*w = 2;return(t_id);}
1698 >
1699 > #define  check_verts(p,v0,v1,v2,o,w,t_id) \
1700 >    if(EQUAL_VEC3(v0,p)){ *o = ON_V; *w = 0; return(t_id); } \
1701 >    if(EQUAL_VEC3(v1,p)){ *o = ON_V; *w = 1; return(t_id); } \
1702 >    if(EQUAL_VEC3(v2,p)){ *o = ON_V; *w = 2; return(t_id);}
1703 >
1704 >
1705 >
1706 > find_neighbor(argptr,qt,f,doneptr)
1707 > int *argptr;
1708 > QUADTREE qt;
1709 > FUNC f;
1710 > int *doneptr;
1711 > {
1712 >  S_ID s_id;
1713 >  int i;
1714 >  S_ID *optr;
1715 >  
1716 >  if(QT_IS_EMPTY(qt))
1717 >    return;
1718 >  else
1719 >    if(QT_IS_TREE(qt))
1720      {
1721 <      if(on2)
1721 >      for(i=0;i < 4; i++)
1722        {
1723 <        *onptr = ON_P;
1724 <        *whichptr = 0;
1723 >        find_neighbor(argptr,QT_NTH_CHILD(qt,i),f,doneptr);
1724 >        if(*doneptr)
1725 >          return;
1726        }
1727 <      else
1727 >    }
1728 >    else
1729 >    {
1730 >      optr = qtqueryset(qt);
1731 >      for(i = QT_SET_CNT(optr); i > 0; i--)
1732        {
1733 <        *onptr = ON_E;
1734 <        *whichptr = 1;
1733 >        s_id = QT_SET_NEXT_ELEM(optr);
1734 >        if(s_id != ((S_ARGS *)argptr)->s_id)
1735 >        {
1736 >          *doneptr = TRUE;
1737 >          ((S_ARGS *)argptr)->n_id = s_id;
1738 >          return;
1739 >        }
1740        }
1175      return(t_id);
1741      }
1742 <   if(on2)
1743 <   {
1744 <     *onptr = ON_E;
1745 <     *whichptr = 2;
1746 <      return(t_id);
1747 <   }
1748 <    *onptr = IN_T;
1749 <    return(t_id);
1742 > }
1743 >
1744 > smInsert_samp(sm,s_id,p,nptr)
1745 > SM *sm;
1746 > S_ID s_id;
1747 > FVECT p;
1748 > S_ID *nptr;
1749 > {
1750 >  FVECT tpt;
1751 >  FUNC f;
1752 >  S_ARGS args;
1753 >
1754 >  F_FUNC(f) = insert_samp;
1755 >  args.s_id = s_id;
1756 >  if(nptr)
1757 >  {
1758 >    args.n_id = 1;
1759 >    F_FUNC2(f) = find_neighbor;
1760    }
1761 <  return(INVALID);
1761 >  else
1762 >    args.n_id = 0;
1763 >  F_ARGS(f) = (int *)(&args);
1764 >  stInsert_samp(SM_LOCATOR(sm),p,f);
1765 >
1766 >  if(nptr)
1767 >    *nptr = args.n_id;
1768   }
1769  
1770 + /* Assumed p is normalized to view sphere */
1771   int
1772 < smPointLocateTri(sm,p,onptr,whichptr)
1772 > smSample_locate_tri(sm,p,s_id,onptr,whichptr,nptr)
1773   SM *sm;
1774   FVECT p;
1775 + S_ID s_id;
1776   int *onptr,*whichptr;
1777 + S_ID *nptr;
1778   {
1195  QUADTREE qt,*optr;
1196  FVECT tpt;
1779    int tri_id;
1780 +  FVECT tpt;
1781 +  TRI *tri;
1782  
1783 <  VSUB(tpt,p,SM_VIEW_CENTER(sm));
1784 <  qt = smPointLocateCell(sm,tpt);
1785 <  if(QT_IS_EMPTY(qt))
1786 <    return(INVALID);  
1787 <
1788 <  optr = qtqueryset(qt);
1789 <  tri_id = smTri_in_set(sm,tpt,optr,onptr,whichptr);
1783 >  if(SM_S_NTH_QT(sm,s_id) == EMPTY)
1784 >    smInsert_samp(sm,s_id,p,nptr);
1785 > #ifdef DEBUG
1786 >  if(*nptr == INVALID)
1787 >    error(CONSISTENCY,"smSample_locate_tri: unable to find sample\n");
1788 > #endif
1789 >  
1790 >  /* Get the triangle adjacent to a sample in qt, or if there are no other
1791 >     samples in qt than the one just inserted, look at siblings-from parent
1792 >     node
1793 >  */
1794  
1795 <  return(tri_id);
1795 >  tri_id = SM_NTH_VERT(sm,*nptr);
1796 > #ifdef DEBUG
1797 >  if(tri_id == INVALID)
1798 >      error(CONSISTENCY,"smSample_locate_tri: Invalid tri_id\n");
1799 > #endif
1800 >
1801 >  tri = SM_NTH_TRI(sm,tri_id);
1802 > #ifdef DEBUG
1803 >  Tcnt++;
1804 >    if(!T_IS_VALID(tri))
1805 >      error(CONSISTENCY,"smSample_locate_tri: Invalid tri\n");
1806 > #endif
1807 >  tri_id = smWalk(sm,p,tri_id,tri,onptr,whichptr,-1,0.0,NULL,NULL);
1808 > #ifdef DEBUG
1809 >    if(tri_id == INVALID)
1810 >      error(CONSISTENCY,"smSample_locate_tri: unable to find tri on walk\n");
1811 > #endif
1812 >    return(tri_id);
1813   }
1814  
1210
1815   /*
1816     Determine whether this sample should:
1817     a) be added to the mesh by subdividing the triangle
# Line 1237 | Line 1841 | int *onptr,*whichptr;
1841       This attempts to throw out points that should be occluded  
1842   */
1843   int
1844 < smTest_sample(sm,tri_id,dir,p,rptr)
1844 > smTest_samp(sm,tri_id,dir,p,rptr,ns,n0,n1,n2,ds,d0,on,which)
1845     SM *sm;
1846     int tri_id;
1847     FVECT dir,p;
1848 <   int *rptr;
1848 >   S_ID *rptr;
1849 >   FVECT ns,n0,n1,n2;
1850 >   double ds,d0;
1851 >   int on,which;
1852   {
1853      TRI *tri;
1854 <    double d,d2,dnear,dspt,d01,d12,d20,s0,s1,s2,ds,dv;
1855 <    int vid[3],i,nearid,norm,dcnt,bcnt;
1856 <    FVECT diff[3],spt,npt;
1854 >    double d,d2,dnear,dspt,d01,d12,d20,s0,s1,s2;
1855 >    double dp,dp1;
1856 >    S_ID vid[3];
1857 >    int i,norm,dcnt,bcnt;
1858 >    FVECT diff[3],spt,npt,n;
1859 >    FVECT nearpt;
1860  
1861      *rptr = INVALID;
1862      bcnt = dcnt = 0;
# Line 1259 | Line 1869 | smTest_sample(sm,tri_id,dir,p,rptr)
1869      {
1870        if(SM_BASE_ID(sm,vid[i]))
1871          bcnt++;
1872 <      if(SM_DIR_ID(sm,vid[i]))
1873 <        dcnt++;
1874 <      
1872 >      else
1873 >        if(SM_DIR_ID(sm,vid[i]))
1874 >          dcnt++;
1875      }
1876 <    /* TEST 1: If the new sample is close in ws, and close in the spherical
1877 <       projection to one of the triangle vertex samples
1878 <    */
1879 < #if 0
1880 <    norm = FALSE;
1881 <      VSUB(spt,p,SM_VIEW_CENTER(sm));
1882 <      ds = DOT(spt,spt);
1883 <      dnear = FHUGE;
1884 <      for(i=0; i<3; i++)
1885 <        {
1886 <          d = DOT(diff[i],diff[i])/ds;
1887 <          if(d < dnear)
1888 <            {
1889 <              dnear = d;
1890 <              nearid=vid[i];
1891 <            }
1892 <        }
1893 <
1894 <      if(dnear <=  SM_MIN_SAMP_DIFF*SM_MIN_SAMP_DIFF)
1895 <        {
1896 <          /* Pick the point with dir closest to that of the canonical view
1897 <             if it is the new sample: mark existing point for deletion
1898 <             */
1899 <          if(SM_BASE_ID(sm,nearid))
1900 <          {
1901 <            *rptr = nearid;
1902 <            return(TRUE);
1903 <          }
1904 <          if(SM_DIR_ID(sm,nearid))
1905 <            return(FALSE);
1906 <          if(!dir)
1297 <          {
1298 <            *rptr = nearid;
1299 <            return(TRUE);
1300 <          }
1301 <          normalize(spt);
1302 <          norm = TRUE;
1303 <          VSUB(npt,SM_NTH_WV(sm,nearid),SM_VIEW_CENTER(sm));
1304 <          normalize(npt);
1305 <          d = fdir2diff(SM_NTH_W_DIR(sm,nearid), npt);
1306 <          d2 = 2. - 2.*DOT(dir,spt);
1307 <          /* The existing sample is a better sample:punt */
1308 <         if(d2 > d)
1309 <            return(FALSE);
1310 <         else
1311 <         {
1312 <        /* The new sample is better: mark the existing one
1313 <                   for deletion after the new one is added*/
1314 <        *rptr = nearid;
1876 >    if(on == ON_P)
1877 >    {
1878 >      FVECT edir;
1879 >     /* Pick the point with dir closest to that of the canonical view
1880 >         if it is the new sample: mark existing point for deletion
1881 >         */
1882 >      if(!dir)
1883 >        return(FALSE);
1884 >      if(SM_BASE_ID(sm,vid[which]))
1885 >        return(FALSE);
1886 >      if(SM_DIR_ID(sm,vid[which]))
1887 >      {
1888 >        *rptr = vid[which];
1889 >        return(TRUE);
1890 >      }
1891 >  decodedir(edir,SM_NTH_W_DIR(sm,vid[which]));
1892 >      if(which == 0)
1893 >        d = DOT(edir,n0);
1894 >      else
1895 >        if(which == 1)
1896 >          d = DOT(edir,n1);
1897 >        else
1898 >          d = DOT(edir,n2);
1899 >      d2 = DOT(dir,ns);
1900 >      /* The existing sample is a better sample:punt */
1901 >      if(d2 >= d)
1902 >        return(FALSE);
1903 >      else
1904 >      {
1905 >        /* The new sample is better: mark existing one for deletion*/
1906 >        *rptr = vid[which];
1907          return(TRUE);
1908 <        }
1909 <     }  
1318 <    
1319 <  /* TEST 3: If the spherical projection of new is < S_REPLACE_EPS
1320 <     from a base point: Edge length is constrained to subtend <45 degrees:
1321 <     original base mesh edges are approx 32 degrees- so have about 13 degrees
1322 <     to work in: S_REPLACE_EPS is the square of the radian value
1323 <  */
1908 >      }
1909 >    }  
1910  
1911 <    if(bcnt)
1912 <    {  
1913 <        dnear = FHUGE;
1914 <        if(!norm)
1915 <          normalize(spt);
1911 > /* Now test if sample epsilon is within circumcircle of enclosing tri
1912 >   if not punt:
1913 >   */
1914 >    {
1915 >      double ab[3],ac[3],r[3];
1916  
1917 <        for(i=0; i<3; i++)
1917 >      VSUB(ab,n2,n1);
1918 >      VSUB(ac,n0,n2);
1919 >      VCROSS(r,ab,ac);
1920 >      dp = DOT(r,n1);
1921 >      dp1 = DOT(r,ns);
1922 >      if(dp1*dp1 + COS_VERT_EPS*DOT(r,r) < dp*dp)
1923          {
1924 <            if(!SM_BASE_ID(sm,vid[i]))
1925 <               continue;
1926 <            VSUB(npt,SM_NTH_WV(sm,vid[i]),SM_VIEW_CENTER(sm));
1927 <            d = DIST_SQ(npt,spt);
1337 <            if(dnear <=  SM_MIN_SAMP_DIFF*SM_MIN_SAMP_DIFF && d< near)
1338 <               {
1339 <                   dnear = d;
1340 <                   nearid = vid[i];
1341 <               }
1924 > #ifdef DEBUG
1925 >          eputs("smTest_samp:rejecting samp,cant guarantee not within eps\n");
1926 > #endif
1927 >          return(FALSE);
1928          }
1343        if(dnear != FHUGE)
1344        {
1345            /* add new and mark the base for deletion */
1346            *rptr = nearid;
1347            return(TRUE);
1348        }
1929      }
1930 < #else
1930 > #ifdef DEBUG
1931 > #if DEBUG > 1
1932 >    if(on == ON_E)
1933      {
1934 <      FVECT nearpt;
1935 <    dnear = FHUGE;
1936 <    VSUB(spt,p,SM_VIEW_CENTER(sm));
1355 <    ds = DOT(spt,spt);
1356 <    normalize(spt);
1934 >      TRI *topp;
1935 >      int opp_id;
1936 >      FVECT vopp;
1937  
1938 <    for(i=0; i<3; i++)
1939 <    {
1940 <        
1941 <      VSUB(npt,SM_NTH_WV(sm,vid[i]),SM_VIEW_CENTER(sm));
1942 <      
1943 <      if(!SM_BASE_ID(sm,vid[i]) && !SM_DIR_ID(sm,vid[i]))
1944 <        normalize(npt);
1945 <
1946 <      d = DIST_SQ(npt,spt);
1947 <      if(d < SM_MIN_SAMP_DIFF*SM_MIN_SAMP_DIFF && d < dnear)
1948 <        {
1949 <          dnear = d;
1950 <          nearid = vid[i];
1951 <          VCOPY(nearpt,npt);
1952 <        }
1953 <
1938 >      topp = SM_NTH_TRI(sm,T_NTH_NBR(tri,which));
1939 >      opp_id = T_NTH_V(topp, T_NTH_NBR_PTR(tri_id,topp));
1940 >      /* first check if valid tri created*/
1941 >      VSUB(vopp,SM_NTH_WV(sm,opp_id),SM_VIEW_CENTER(sm));
1942 >      normalize(vopp);
1943 >      VCROSS(n,ns,vopp);
1944 >      normalize(n);
1945 >      if(which==2)
1946 >        dp = DOT(n,n0) * DOT(n,n1);
1947 >      else
1948 >        if(which==0)
1949 >          dp = DOT(n,n1) * DOT(n,n2);
1950 >        else
1951 >          dp = DOT(n,n2) * DOT(n,n0);
1952 >      if(dp > 0.0)
1953 >      {
1954 >        eputs("smInsert_samp_mesh: couldn't split edge: rejecting\n");
1955 >        /* test epsilon for edge */
1956 >        return(FALSE);
1957 >      }
1958      }
1959 <    if(dnear != FHUGE)
1960 <    {
1377 <          /* Pick the point with dir closest to that of the canonical view
1378 <             if it is the new sample: mark existing point for deletion
1379 <             */
1380 <          if(SM_BASE_ID(sm,nearid))
1381 <          {
1382 <            *rptr = nearid;
1383 <            return(TRUE);
1384 <          }
1385 <          if(SM_DIR_ID(sm,nearid))
1386 <            return(FALSE);
1387 <          if(!dir)
1388 <          {
1389 <            *rptr = nearid;
1390 <            return(TRUE);
1391 <          }
1392 <          d = fdir2diff(SM_NTH_W_DIR(sm,nearid), nearpt);
1393 <          d2 = 2. - 2.*DOT(dir,spt);
1394 <          /* The existing sample is a better sample:punt */
1395 <         if(d2 > d)
1396 <            return(FALSE);
1397 <         else
1398 <         {
1399 <        /* The new sample is better: mark the existing one
1400 <                   for deletion after the new one is added*/
1401 <        *rptr = nearid;
1402 <        return(TRUE);
1403 <        }
1404 <    }  
1405 <    }
1406 < #endif  
1959 > #endif
1960 > #endif
1961    /* TEST 4:
1962       If the addition of the new sample point would introduce a "puncture"
1963       or cause new triangles with large depth differences:dont add:    
# Line 1411 | Line 1965 | smTest_sample(sm,tri_id,dir,p,rptr)
1965      if(bcnt || dcnt)
1966         return(TRUE);
1967  
1968 <    /* If the new point is in front of the existing points- add */
1415 <    dv = DIST_SQ(SM_NTH_WV(sm,vid[0]),SM_VIEW_CENTER(sm));
1416 <    if(ds < dv)
1968 >    if(ds < d0)
1969        return(TRUE);
1970  
1971      d01 = DIST_SQ(SM_NTH_WV(sm,vid[1]),SM_NTH_WV(sm,vid[0]));
# Line 1441 | Line 1993 | smTest_sample(sm,tri_id,dir,p,rptr)
1993    /* TEST 5:
1994       Check to see if triangle is relatively large and should therefore
1995       be subdivided anyway.
1996 <     */
1997 <    dv *= DIST_SQ(SM_NTH_WV(sm,vid[1]),SM_VIEW_CENTER(sm));
1998 <    dv *= DIST_SQ(SM_NTH_WV(sm,vid[2]),SM_VIEW_CENTER(sm));
1999 <    if (d01*d12*d20/dv > S_REPLACE_TRI)
1996 >     */  
1997 >    d0 = d0*d0*DIST_SQ(SM_NTH_WV(sm,vid[1]),SM_VIEW_CENTER(sm));
1998 >    d0 *= DIST_SQ(SM_NTH_WV(sm,vid[2]),SM_VIEW_CENTER(sm));
1999 >    if (d01*d12*d20/d0 > S_REPLACE_TRI)
2000          return(TRUE);
2001              
2002      return(FALSE);
2003   }
2004 <
2005 < smReplace_base_samp(sm,b_id,s_id,tri,t_id,which)
1454 < SM *sm;
1455 < int b_id,s_id;
1456 < TRI *tri;
1457 < int t_id,which;
1458 < {
1459 <  TRI *t;
1460 <
1461 <  SM_NTH_VERT(sm,s_id) = t_id;
1462 <  T_NTH_V(tri,which) = s_id;
1463 <  if(!(SM_BASE_ID(sm,T_NTH_V(tri,(which+1)%3)) ||
1464 <       SM_BASE_ID(sm,T_NTH_V(tri,(which+2)%3))))
1465 <    SM_CLR_NTH_T_BASE(sm,t_id);
1466 <  if(!SM_IS_NTH_T_NEW(sm,t_id))
1467 <  {
1468 <    smNew_tri_cnt++;
1469 <    SM_SET_NTH_T_NEW(sm,t_id);
1470 <  }
1471 <  t_id = smTri_next_ccw_nbr(sm,tri,b_id);
1472 <  while(t_id != INVALID)
1473 <  {
1474 <    t = SM_NTH_TRI(sm,t_id);
1475 <    which = T_WHICH_V(t,b_id);
1476 <    T_NTH_V(t,which) = s_id;
1477 <    /* Check if still a base triangle */
1478 <    if(!(SM_BASE_ID(sm,T_NTH_V(t,(which+1)%3)) ||
1479 <         SM_BASE_ID(sm,T_NTH_V(t,(which+2)%3))))
1480 <      SM_CLR_NTH_T_BASE(sm,t_id);
1481 <     if(!SM_IS_NTH_T_NEW(sm,t_id))
1482 <        {
1483 <          smNew_tri_cnt++;
1484 <          SM_SET_NTH_T_NEW(sm,t_id);
1485 <        }
1486 <     t_id = smTri_next_ccw_nbr(sm,t,b_id);
1487 <    }
1488 < }
1489 <
1490 < int
1491 < smReplace_samp(sm,c,dir,p,s_id,t_id,o_id,on,which)
2004 > S_ID
2005 > smReplace_samp(sm,c,dir,p,np,s_id,t_id,o_id,on,which)
2006       SM *sm;
2007       COLR c;
2008 <     FVECT dir,p;
2009 <     int s_id,t_id,o_id,on,which;
2008 >     FVECT dir,p,np;
2009 >     S_ID s_id;
2010 >     int t_id;
2011 >     S_ID o_id;
2012 >     int on,which;
2013   {
2014 <  int v_id,tri_id;
2014 >  S_ID v_id;
2015 >  int tri_id;
2016    TRI *t,*tri;
2017  
2018    tri = SM_NTH_TRI(sm,t_id);
# Line 1502 | Line 2020 | smReplace_samp(sm,c,dir,p,s_id,t_id,o_id,on,which)
2020  
2021    /* If it is a base id, need new sample */
2022    if(SM_BASE_ID(sm,v_id))
2023 <  {
2024 <    sInit_samp(SM_SAMP(sm),s_id,c,dir,p,o_id);
2025 <    smReplace_base_samp(sm,v_id,s_id,tri,t_id,which);
2026 <    return(s_id);
2027 <  }
2023 >    return(v_id);
2024 >
2025 >  if(EQUAL_VEC3(p,SM_NTH_WV(sm,v_id)))
2026 >    return(INVALID);
2027 >
2028    if(dir)
2029 <  {
2030 <    /* If world point */
2031 <    /* if existing point is a dir: leave */
2032 <    if(SM_DIR_ID(sm,v_id))
2033 <      return(INVALID);
2034 <    if(on == ON_V)
2029 >  {    /* If world point */
2030 >     FVECT odir,no;
2031 >     double d,d1;
2032 >     int dir_id;
2033 >    /* if existing point is a not a dir, compare directions */
2034 >    if(!(dir_id = SM_DIR_ID(sm,v_id)))
2035      {
2036 +      decodedir(odir, SM_NTH_W_DIR(sm,v_id));  
2037 +      VSUB(no,SM_NTH_WV(sm,v_id),SM_VIEW_CENTER(sm));
2038 +      d = DOT(dir,np);
2039 +      d1 = DOT(odir,np);
2040 +     /* The existing sample is a better sample:punt */
2041 +    }
2042 +    if(dir_id || (d*d*DOT(no,no) > d1*d1))
2043 +    {
2044 +      /* The new sample has better information- replace values */
2045        sInit_samp(SM_SAMP(sm),v_id,c,dir,p,o_id);
2046 <
2047 <     if(!SM_IS_NTH_T_NEW(sm,t_id))
1521 <        {
1522 <          smNew_tri_cnt++;
1523 <          SM_SET_NTH_T_NEW(sm,t_id);
1524 <        }
2046 >      if(!SM_IS_NTH_T_NEW(sm,t_id))
2047 >         SM_SET_NTH_T_NEW(sm,t_id);
2048        tri_id = smTri_next_ccw_nbr(sm,tri,v_id);
2049        while(tri_id != t_id)
2050        {
2051          t = SM_NTH_TRI(sm,tri_id);
2052 <     if(!SM_IS_NTH_T_NEW(sm,tri_id))
1530 <        {
1531 <          smNew_tri_cnt++;
2052 >        if(!SM_IS_NTH_T_NEW(sm,tri_id))
2053            SM_SET_NTH_T_NEW(sm,tri_id);
1533        }
1534
2054          tri_id = smTri_next_ccw_nbr(sm,t,v_id);
2055        }
1537      return(v_id);
2056      }
1539    /* on == ON_P */
1540    else
1541   {
1542     FVECT spt,npt;
1543     double d,d2;
1544
1545    /* Need to check which sample is preferable */
1546     VSUB(spt,p,SM_VIEW_CENTER(sm));
1547     normalize(spt);
1548        
1549     VSUB(npt,SM_NTH_WV(sm,v_id),SM_VIEW_CENTER(sm));
1550     normalize(npt);
1551     d = fdir2diff(SM_NTH_W_DIR(sm,v_id), npt);
1552     d2 = 2. - 2.*DOT(dir,spt);
1553      /* The existing sample is a better sample:punt */
1554     if(d2 < d)
1555       {
1556         /* The new sample has better information- replace values */
1557         sInit_samp(SM_SAMP(sm),v_id,c,dir,p,o_id);
1558         if(!SM_IS_NTH_T_NEW(sm,t_id))
1559            {
1560              smNew_tri_cnt++;
1561              SM_SET_NTH_T_NEW(sm,t_id);
1562            }
1563
1564         tri_id = smTri_next_ccw_nbr(sm,tri,v_id);
1565         while(tri_id != t_id)
1566           {
1567             t = SM_NTH_TRI(sm,tri_id);
1568             if(!SM_IS_NTH_T_NEW(sm,tri_id))
1569            {
1570              smNew_tri_cnt++;
1571              SM_SET_NTH_T_NEW(sm,tri_id);
1572            }
1573             tri_id = smTri_next_ccw_nbr(sm,t,v_id);
1574           }
1575       }
2057      return(v_id);
1577   }
2058    }
2059 <  else
2060 <    { /* New point is a dir */
2061 <      return(INVALID);
1582 <    }
2059 >  else  /* New point is a dir */
2060 >  return(INVALID);
2061 >
2062   }
2063  
2064 < int
2064 > S_ID
2065   smAlloc_samp(sm)
2066   SM *sm;
2067   {
2068 <  int s_id,replaced,cnt;
2068 >  S_ID s_id;
2069 >  int replaced,cnt;
2070    SAMP *s;
2071    FVECT p;
2072  
2073    s = SM_SAMP(sm);
2074    s_id = sAlloc_samp(s,&replaced);
1595
1596  cnt=0;
2075    while(replaced)
2076    {
2077 <    if(smRemoveVertex(sm,s_id))
1600 <      break;
2077 >    smRemoveVertex(sm,s_id);
2078      s_id = sAlloc_samp(s,&replaced);
1602    cnt++;
1603    if(cnt > S_MAX_SAMP(s))
1604      error(CONSISTENCY,"smAlloc_samp():unable to find free samp\n");
2079    }
2080    return(s_id);
2081   }
# Line 1620 | Line 2094 | SM *sm;
2094          b) coincide with existing edge
2095          c) Fall in existing triangle
2096   */
2097 < int
2097 > S_ID
2098   smAdd_samp(sm,c,dir,p,o_id)
2099     SM *sm;
2100     COLR c;
2101     FVECT dir,p;
2102 <   int o_id;
2102 >   S_ID o_id;
2103   {
2104 <  int t_id,s_id,r_id,on,which,n_id;
2105 <  double d;
2106 <  FVECT wpt;
2104 >  int t_id,on,which,n_id;
2105 >  S_ID s_id,r_id,nbr_id;
2106 >  double ds,d0;
2107 >  FVECT wpt,ns,n0,n1,n2;
2108 >  QUADTREE qt,parent;
2109 >  TRI *t;
2110  
2111    r_id = INVALID;
2112 <
2112 >  nbr_id = INVALID;
2113    /* Must do this first-as may change mesh */
2114    s_id = smAlloc_samp(sm);
2115 +  SM_S_NTH_QT(sm,s_id)= EMPTY;
2116    /* If sample is a world space point */
2117    if(p)
2118    {
2119 +    VSUB(ns,p,SM_VIEW_CENTER(sm));
2120 +    ds = normalize(ns);
2121      while(1)
2122      {
2123 <      t_id = smPointLocateTri(sm,p,&on,&which);
2123 >      t_id = smSample_locate_tri(sm,ns,s_id,&on,&which,&nbr_id);
2124 >      qt = SM_S_NTH_QT(sm,s_id);
2125        if(t_id == INVALID)
2126        {
2127   #ifdef DEBUG
2128            eputs("smAddSamp(): unable to locate tri containing sample \n");
2129   #endif
2130            smUnalloc_samp(sm,s_id);
2131 +          qtremovelast(qt,s_id);
2132            return(INVALID);
2133          }
2134        /* If spherical projection coincides with existing sample: replace */
2135 <      if((on == ON_V || on == ON_P))
2135 >      if(on == ON_V)
2136        {
2137 <          n_id = smReplace_samp(sm,c,dir,p,s_id,t_id,o_id,on,which);
2138 <          if(n_id!= s_id)
2139 <             smUnalloc_samp(sm,s_id);
2137 >          n_id = smReplace_samp(sm,c,dir,p,ns,s_id,t_id,o_id,on,which);
2138 >          if(n_id != s_id)
2139 >          {
2140 >            smUnalloc_samp(sm,s_id);
2141 > #if 0
2142 >            fprintf(stderr,"Unallocating sample %d\n",s_id);
2143 > #endif
2144 >            qtremovelast(qt,s_id);
2145 >          }
2146            return(n_id);
2147        }
2148 <      if((!(smTest_sample(sm,t_id,dir,p,&r_id))))
2149 <     {
2148 >      t = SM_NTH_TRI(sm,t_id);
2149 >      VSUB(n0,SM_NTH_WV(sm,T_NTH_V(t,0)),SM_VIEW_CENTER(sm));
2150 >      d0 = normalize(n0);
2151 >      VSUB(n1,SM_NTH_WV(sm,T_NTH_V(t,1)),SM_VIEW_CENTER(sm));
2152 >      normalize(n1);
2153 >      VSUB(n2,SM_NTH_WV(sm,T_NTH_V(t,2)),SM_VIEW_CENTER(sm));
2154 >      normalize(n2);
2155 >      if((!(smTest_samp(sm,t_id,dir,p,&r_id,ns,n0,n1,n2,ds,d0,on,which))))
2156 >      {
2157          smUnalloc_samp(sm,s_id);
2158 + #if 0
2159 +            fprintf(stderr,"Unallocating sample %d\n",s_id);
2160 + #endif
2161 +        qtremovelast(qt,s_id);
2162          return(INVALID);
2163        }
2164        if(r_id != INVALID)
2165        {
2166          smRemoveVertex(sm,r_id);
2167 <        sDelete_samp(SM_SAMP(sm),r_id);
2167 >        if(nbr_id == r_id)
2168 >        {
2169 >          qtremovelast(qt,s_id);
2170 >          SM_S_NTH_QT(sm,s_id) = EMPTY;
2171 >        }
2172        }
2173        else
2174          break;
# Line 1673 | Line 2176 | smAdd_samp(sm,c,dir,p,o_id)
2176      /* If sample is being added in the middle of the sample array: tone
2177         map individually
2178         */
2179 +    /* If not base or sky point, add distance from center to average*/  
2180 +    smDist_sum += 1.0/ds;
2181      /* Initialize sample */
2182      sInit_samp(SM_SAMP(sm),s_id,c,dir,p,o_id);
2183 <    
2183 >    smInsert_samp_mesh(sm,s_id,t_id,ns,n0,n1,n2,on,which);
2184    }
2185      /* If sample is a direction vector */
2186      else
2187      {
2188        VADD(wpt,dir,SM_VIEW_CENTER(sm));
2189 +      /* Allocate space for a sample and initialize */
2190        while(1)
2191 <        {
2192 <            t_id = smPointLocateTri(sm,wpt,&on,&which);
2193 <            if(t_id == INVALID)
1688 <               {
2191 >      {
2192 >        t_id = smSample_locate_tri(sm,dir,s_id,&on,&which,&nbr_id);
2193 >        qt = SM_S_NTH_QT(sm,s_id);
2194   #ifdef DEBUG
2195 <                   eputs("smAddSamp(): can'tlocate tri containing samp\n");
2195 >        if(t_id == INVALID)
2196 >        {
2197 >          eputs("smAddSamp():can't locate tri containing samp\n");
2198 >          smUnalloc_samp(sm,s_id);
2199 >          qtremovelast(qt,s_id);
2200 >          return(INVALID);
2201 >        }
2202   #endif
2203 <                   smUnalloc_samp(sm,s_id);
2204 <                   return(INVALID);
2205 <               }
2206 <            if(on == ON_V || on == ON_P)
2207 <            {
2208 <                n_id=smReplace_samp(sm,c,NULL,wpt,s_id,t_id,o_id,on,which);
2209 <                if(n_id !=s_id)
2210 <                   smUnalloc_samp(sm,s_id);
2211 <                return(n_id);
2212 <            }
2213 <            if((!(smTest_sample(sm,t_id,NULL,wpt,&r_id))))
2203 >        if(on == ON_V)
2204 >        {
2205 >          smUnalloc_samp(sm,s_id);
2206 >          qtremovelast(qt,s_id);
2207 >          return(INVALID);
2208 >        }
2209 >        t = SM_NTH_TRI(sm,t_id);
2210 >        VSUB(n0,SM_NTH_WV(sm,T_NTH_V(t,0)),SM_VIEW_CENTER(sm));
2211 >        d0 = normalize(n0);
2212 >        VSUB(n1,SM_NTH_WV(sm,T_NTH_V(t,1)),SM_VIEW_CENTER(sm));
2213 >        normalize(n1);
2214 >        VSUB(n2,SM_NTH_WV(sm,T_NTH_V(t,2)),SM_VIEW_CENTER(sm));
2215 >        normalize(n2);
2216 >        if(!smTest_samp(sm,t_id,NULL,wpt,&r_id,dir,n0,n1,n2,1.0,d0,on,which))
2217 >        {
2218 >          smUnalloc_samp(sm,s_id);
2219 >          qtremovelast(qt,s_id);
2220 >          return(INVALID);
2221 >        }
2222 >        if(r_id != INVALID)
2223 >        {
2224 >          smRemoveVertex(sm,r_id);
2225 >          if(nbr_id == r_id)
2226              {
2227 <                smUnalloc_samp(sm,s_id);
2228 <                return(INVALID);
2227 >              qtremovelast(qt,s_id);
2228 >              SM_S_NTH_QT(sm,s_id) = EMPTY;
2229              }
1707            if(r_id != INVALID)
1708            {
1709                smRemoveVertex(sm,r_id);
1710                sDelete_samp(SM_SAMP(sm),r_id);
1711            }
1712            else
1713               break;
2230          }
2231 <      /* Allocate space for a sample and initialize */
2231 >        else
2232 >          break;
2233 >      }
2234        sInit_samp(SM_SAMP(sm),s_id,c,NULL,wpt,o_id);
2235 +      smInsert_samp_mesh(sm,s_id,t_id,dir,n0,n1,n2,on,which);
2236      }
1718  if(!SM_DIR_ID(sm,s_id))
1719    {
1720      /* If not a base or sky point, add distance from the
1721         viewcenter to average*/
1722      d = DIST(SM_NTH_WV(sm,s_id),SM_VIEW_CENTER(sm));
1723      smDist_sum += 1.0/d;
1724    }
1725    smInsert_samp(sm,s_id,t_id,on,which);
1726
2237      return(s_id);
2238   }
2239  
# Line 1738 | Line 2248 | smAdd_samp(sm,c,dir,p,o_id)
2248   * New sample representation will be output in next call to smUpdate().
2249   * If the point is a sky point: then v=NULL
2250   */
2251 + #define FVECT_TO_SFLOAT(p) \
2252 +        (p[0]=(SFLOAT)p[0],p[1]=(SFLOAT)p[1],p[2]=(SFLOAT)p[2])
2253   int
2254   smNewSamp(c,dir,p)
2255   COLR c;
2256   FVECT dir;
2257   FVECT p;
2258   {
2259 <    int s_id;
2259 >    S_ID s_id;
2260  
2261      /* First check if this the first sample: if so initialize mesh */
2262      if(SM_NUM_SAMP(smMesh) == 0)
# Line 1752 | Line 2264 | FVECT p;
2264        smInit_sm(smMesh,odev.v.vp);
2265        sClear_all_flags(SM_SAMP(smMesh));
2266      }
2267 +    FVECT_TO_SFLOAT(p);
2268 +    FVECT_TO_SFLOAT(dir);
2269      /* Add the sample to the mesh */
2270      s_id = smAdd_samp(smMesh,c,dir,p,INVALID);
2271  
2272 <    return(s_id);
2272 >    return((int)s_id);
2273      
2274   }    
2275 < int
2275 > S_ID
2276   smAdd_base_vertex(sm,v)
2277     SM *sm;
2278     FVECT v;
2279   {
2280 <  int id;
2280 >  S_ID id;
2281  
2282    /* First add coordinate to the sample array */
2283    id = sAdd_base_point(SM_SAMP(sm),v);
# Line 1787 | Line 2301 | SM *sm;
2301   int type;
2302   {
2303    int i,s_id,tri_id,nbr_id;
2304 <  int p[SM_EXTRA_POINTS],ids[SM_BASE_TRIS];
2305 <  int v0_id,v1_id,v2_id;
2304 >  S_ID p[SM_BASE_POINTS];
2305 >  int ids[SM_BASE_TRIS];
2306 >  S_ID v0_id,v1_id,v2_id;
2307    FVECT d,pt,cntr,v0,v1,v2;
2308    TRI *t;
2309  
2310    /* First insert the base vertices into the sample point array */
2311 <
1797 <  for(i=0; i < SM_EXTRA_POINTS; i++)
2311 >  for(i=0; i < SM_BASE_POINTS; i++)
2312    {
2313      VADD(cntr,icosa_verts[i],SM_VIEW_CENTER(sm));
2314      p[i]  = smAdd_base_vertex(sm,cntr);
2315 +    smInsert_samp(sm,p[i],icosa_verts[i],NULL);
2316    }
2317    /* Create the base triangles */
2318    for(i=0;i < SM_BASE_TRIS; i++)
# Line 1805 | Line 2320 | int type;
2320      v0_id = p[icosa_tris[i][0]];
2321      v1_id = p[icosa_tris[i][1]];
2322      v2_id = p[icosa_tris[i][2]];
2323 <    ids[i] = smAdd_tri(sm, v0_id,v1_id,v2_id);
2323 >    ids[i] = smAdd_tri(sm, v0_id,v1_id,v2_id,&t);
2324      /* Set neighbors */
2325      for(nbr_id=0; nbr_id < 3; nbr_id++)
2326        T_NTH_NBR(SM_NTH_TRI(sm,ids[i]),nbr_id) = icosa_nbrs[i][nbr_id];
2327  
2328    }
1814  for(i=0;i < SM_BASE_TRIS; i++)
1815  {
1816    t = SM_NTH_TRI(sm,ids[i]);
1817    smLocator_add_tri(sm,ids[i],T_NTH_V(t,0),T_NTH_V(t,1),T_NTH_V(t,2));
1818  }
2329   return(1);
2330  
2331   }
2332  
2333 < smRebuild_mesh(sm,v)
2333 > smRebuild(sm,v)
2334     SM *sm;
2335     VIEW *v;
2336   {
2337 <    int i,j,cnt;
2337 >    S_ID s_id;
2338 >    int j,cnt;
2339      FVECT p,ov,dir;
2340      double d;
2341  
2342   #ifdef DEBUG
2343 <    fprintf(stderr,"smRebuild_mesh(): rebuilding....");
2343 >    fprintf(stderr,"smRebuild(): rebuilding....");
2344   #endif
2345      smCull(sm,v,SM_ALL_LEVELS);
2346      /* Clear the mesh- and rebuild using the current sample array */
# Line 1838 | Line 2349 | smRebuild_mesh(sm,v)
2349      /* Calculate the difference between the current and new viewpoint*/
2350      /* Will need to subtract this off of sky points */
2351      VCOPY(ov,SM_VIEW_CENTER(sm));
1841    /* Initialize the mesh to 0 samples and the base triangles */
2352  
2353 +    /* Initialize the mesh to 0 samples and the base triangles */
2354 +    smInit_sm(sm,v->vp);
2355      /* Go through all samples and add in if in the new view frustum, and
2356         the dir is <= 30 degrees from new view
2357       */
2358 <    smInit_sm(sm,v->vp);
1847 <    for(i=0; i < cnt; i++)
2358 >    for(s_id=0; s_id < cnt; s_id++)
2359      {
2360        /* First check if sample visible(conservative approx: if one of tris
2361           attached to sample is in frustum)       */
2362 <      if(!S_IS_FLAG(i))
2362 >      if(!S_IS_FLAG(s_id))
2363          continue;
2364 <
2365 <      /* Next test if direction > 30 from current direction */
2366 <        if(SM_NTH_W_DIR(sm,i)!=-1)
2364 >      
2365 >      /* Next test if direction > 30 degrees from current direction */
2366 >        if(SM_NTH_W_DIR(sm,s_id)!=-1)
2367          {
2368 <            VSUB(p,SM_NTH_WV(sm,i),v->vp);
2369 <            normalize(p);
2370 <            decodedir(dir,SM_NTH_W_DIR(sm,i));
2371 <            d = 2. - 2.*DOT(dir,p);
1861 <            if (d > MAXDIFF2)
2368 >            VSUB(p,SM_NTH_WV(sm,s_id),v->vp);
2369 >            decodedir(dir,SM_NTH_W_DIR(sm,s_id));
2370 >            d = DOT(dir,p);
2371 >            if (d*d < MAXDIFF2*DOT(p,p))
2372                continue;
2373 <            VCOPY(p,SM_NTH_WV(sm,i));
2374 <            smAdd_samp(sm,NULL,dir,p,i);
2373 >            VCOPY(p,SM_NTH_WV(sm,s_id));
2374 >            smAdd_samp(sm,NULL,dir,p,s_id);
2375          }
2376          else
2377          {
2378 <          VSUB(dir,SM_NTH_WV(sm,i),ov);
2379 <          VADD(SM_NTH_WV(sm,i),dir,SM_VIEW_CENTER(sm));
2380 <          smAdd_samp(sm,NULL,dir,NULL,i);
2378 >          /* If the direction is > 45 degrees from current view direction:
2379 >             throw out
2380 >           */
2381 >          VSUB(dir,SM_NTH_WV(sm,s_id),ov);
2382 >          if(DOT(dir,v->vdir) < MAXDIR)
2383 >            continue;
2384 >
2385 >          VADD(SM_NTH_WV(sm,s_id),dir,SM_VIEW_CENTER(sm));
2386 >          smAdd_samp(sm,NULL,dir,NULL,s_id);
2387          }
2388  
2389      }
2390   #ifdef DEBUG
2391 <    fprintf(stderr,"smRebuild_mesh():done\n");
2391 >    fprintf(stderr,"smRebuild():done\n");
2392   #endif
2393   }
2394  
1879 int
1880 intersect_tri_set(t_set,orig,dir,pt)
1881   OBJECT *t_set;
1882   FVECT orig,dir,pt;
1883 {
1884    OBJECT *optr;
1885    int i,t_id,id,base;
1886    int pid0,pid1,pid2;
1887    FVECT p0,p1,p2,p;
1888    TRI *t;
1889    double d,d1;
1890    
1891    optr = QT_SET_PTR(t_set);
1892    for(i = QT_SET_CNT(t_set); i > 0; i--)
1893    {
1894        t_id = QT_SET_NEXT_ELEM(optr);
1895        t = SM_NTH_TRI(smMesh,t_id);
1896        if(!T_IS_VALID(t) || SM_IS_NTH_T_BASE(smMesh,t_id)||
1897               SM_IS_NTH_T_BG(smMesh,t_id))
1898          continue;
1899        pid0 = T_NTH_V(t,0);
1900        pid1 = T_NTH_V(t,1);
1901        pid2 = T_NTH_V(t,2);
1902        VCOPY(p0,SM_NTH_WV(smMesh,pid0));
1903        VCOPY(p1,SM_NTH_WV(smMesh,pid1));
1904        VCOPY(p2,SM_NTH_WV(smMesh,pid2));
1905        if(ray_intersect_tri(orig,dir,p0,p1,p2,p))
1906        {
1907          d =  DIST_SQ(p,p0);
1908          d1 = DIST_SQ(p,p1);
1909          if(d < d1)
1910           {
1911             d1 = DIST_SQ(p,p2);
1912             id = (d1 < d)?pid2:pid0;
1913           }
1914          else
1915            {
1916              d = DIST_SQ(p,p2);
1917              id = (d < d1)? pid2:pid1;
1918            }
1919          if(pt)
1920            VCOPY(pt,p);
1921 #ifdef TEST_DRIVER
1922          Pick_tri = t_id;
1923          Pick_samp = id;
1924          VCOPY(Pick_point[0],p);
1925 #endif
1926          return(id);
1927        }
1928    }
1929    return(-1);
1930 }
1931
2395   /* OS is constrained to be <= QT_MAXCSET : if the set exceeds this, the
2396   results of check_set are conservative
2397   */
2398   int
2399   compare_ids(id1,id2)
2400 < OBJECT *id1,*id2;
2400 > S_ID *id1,*id2;
2401   {
2402    int d;
2403  
# Line 1948 | Line 2411 | OBJECT *id1,*id2;
2411    return(0);
2412   }
2413  
1951 ray_trace_check_set(qt,argptr,fptr)
1952   QUADTREE qt;
1953   RT_ARGS *argptr;
1954   int *fptr;
1955 {
1956    OBJECT tset[QT_MAXSET+1],*tptr;    
1957    double dt,t;
1958    int found;
1959  if(QT_IS_EMPTY(qt))
1960    return;
1961  if(QT_LEAF_IS_FLAG(qt))
1962  {
1963    QT_FLAG_SET_DONE(*fptr);
1964 #if DEBUG
1965       eputs("ray_trace_check_set():Already visited this node:aborting\n");
1966 #endif
1967       return;
1968  }
1969  else
1970    QT_LEAF_SET_FLAG(qt);
2414  
1972  tptr = qtqueryset(qt);
1973  if(QT_SET_CNT(tptr) > QT_MAXSET)
1974    tptr = (OBJECT *)malloc((QT_SET_CNT(tptr)+1)*sizeof(OBJECT));
1975  else
1976    tptr = tset;
1977  if(!tptr)
1978    goto memerr;
1979    
1980  qtgetset(tptr,qt);
1981  /* Must sort */
1982  qsort((void *)(&(tptr[1])),tptr[0],sizeof(OBJECT),compare_ids);
1983  /* Check triangles in set against those seen so far(os):only
1984     check new triangles for intersection (t_set')  
1985     */
1986  check_set_large(tptr,argptr->os);
1987
1988  if(!QT_SET_CNT(tptr))
1989    return;
1990  found = intersect_tri_set(tptr,argptr->orig,argptr->dir,NULL);
1991  if(tptr != tset)
1992    free(tptr);
1993  if(found != INVALID)
1994  {
1995    argptr->t_id = found;
1996    QT_FLAG_SET_DONE(*fptr);
1997  }
1998  return;
1999 memerr:
2000    error(SYSTEM,"ray_trace_check_set():Unable to allocate memory");
2001 }
2002
2003
2004 /*
2005 * int
2006 * smFindSamp(FVECT orig, FVECT dir)
2007 *
2008 * Find the closest sample to the given ray.  Returns sample id, -1 on failure.
2009 * "dir" is assumed to be normalized
2010 */
2011
2012 int
2013 smFindSamp(orig,dir)
2014 FVECT orig,dir;
2015 {
2016  FVECT b,p,o;
2017  OBJECT *ts;
2018  QUADTREE qt;
2019  int s_id,test;
2020  double d;
2021
2022 /*  r is the normalized vector from the view center to the current
2023  *  ray point ( starting with "orig"). Find the cell that r falls in,
2024  *  and test the ray against all triangles stored in the cell. If
2025  *  the test fails, trace the projection of the ray across to the
2026  *  next cell it intersects: iterate until either an intersection
2027  *  is found, or the projection ray is // to the direction. The sample
2028  *  corresponding to the triangle vertex closest to the intersection
2029  *  point is returned.
2030  */
2031  
2032  /* First test if "orig" coincides with the View_center or if "dir" is
2033     parallel to r formed by projecting "orig" on the sphere. In
2034     either case, do a single test against the cell containing the
2035     intersection of "dir" and the sphere
2036   */
2037  /* orig will be updated-so preserve original value */
2038  if(!smMesh)
2039     return;
2040 #ifdef TEST_DRIVER
2041  Picking= TRUE;
2042 #endif  
2043
2044  if(EQUAL_VEC3(orig,SM_VIEW_CENTER(smMesh)))
2045  {
2046    qt = smPointLocateCell(smMesh,dir);
2047    if(QT_IS_EMPTY(qt))
2048      goto Lerror;
2049    ts = qtqueryset(qt);
2050    s_id = intersect_tri_set(ts,orig,dir,p);
2051 #ifdef DEBUG_TEST_DRIVER
2052    VCOPY(Pick_point[0],p);
2053 #endif
2054    return(s_id);
2055  }
2056  d = point_on_sphere(b,orig,SM_VIEW_CENTER(smMesh));
2057  if(EQUAL_VEC3(b,dir))
2058  {
2059    qt = smPointLocateCell(smMesh,dir);
2060    if(QT_IS_EMPTY(qt))
2061      goto Lerror;
2062    ts = qtqueryset(qt);
2063    s_id = intersect_tri_set(ts,orig,dir,p);
2064 #ifdef DEBUG_TEST_DRIVER
2065        VCOPY(Pick_point[0],p);
2066 #endif
2067        return(s_id);
2068  }
2069  if(OPP_EQUAL_VEC3(b,dir))
2070  {
2071    qt = smPointLocateCell(smMesh,orig);
2072    if(QT_IS_EMPTY(qt))
2073      goto Lerror;
2074    ts = qtqueryset(qt);
2075    s_id = intersect_tri_set(ts,orig,dir,p);
2076 #ifdef DEBUG_TEST_DRIVER
2077        VCOPY(Pick_point[0],p);
2078 #endif
2079        if(s_id != INVALID)
2080        {
2081 #ifdef DEBUG
2082    fprintf(stderr,"Front pick returning %d\n",s_id);
2083 #endif
2084          return(s_id);
2085        }
2086        qt = smPointLocateCell(smMesh,dir);
2087        if(QT_IS_EMPTY(qt))
2088          goto Lerror;
2089        ts = qtqueryset(qt);
2090        s_id = intersect_tri_set(ts,orig,dir,p);
2091 #ifdef DEBUG_TEST_DRIVER
2092        VCOPY(Pick_point[0],p);
2093 #endif
2094 #ifdef DEBUG
2095    fprintf(stderr,"Back pick returning %d\n",s_id);
2096 #endif
2097        return(s_id);
2098  }
2099  {
2100    OBJECT t_set[QT_MAXCSET + 1];
2101    RT_ARGS rt;
2102    FUNC func;
2103
2104    /* Test each of the root triangles against point id */
2105    QT_CLEAR_SET(t_set);
2106    VSUB(o,orig,SM_VIEW_CENTER(smMesh));
2107    ST_CLEAR_FLAGS(SM_LOCATOR(smMesh));
2108    rt.t_id = -1;
2109    rt.os = t_set;
2110    VCOPY(rt.orig,orig);
2111    VCOPY(rt.dir,dir);
2112    F_FUNC(func) = ray_trace_check_set;
2113    F_ARGS(func) = (int *)(&rt);
2114    stTrace_ray(SM_LOCATOR(smMesh),o,dir,func);
2115    s_id = rt.t_id;
2116
2117  }    
2118  return(s_id);
2119
2120 Lerror:
2121 #ifdef DEBUG
2122  eputs("smFindSamp(): point not found");
2123 #endif  
2124  return(INVALID);
2125
2126 }
2127
2415   null_func(argptr,root,qt,n)
2416       int *argptr;
2417       int root;
# Line 2134 | Line 2421 | null_func(argptr,root,qt,n)
2421    /* do nothing */
2422   }
2423  
2424 < mark_active_tris(argptr,root,qt,n)
2424 > mark_active_samples(argptr,root,qt,n)
2425       int *argptr;
2426       int root;
2427       QUADTREE qt;
2428       int n;
2429   {
2430 <  OBJECT *os,*optr;
2431 <  register int i,t_id;
2430 >  S_ID *os,s_id;
2431 >  register int i,t_id,tri_id;
2432    TRI *tri;
2433  
2434    if(QT_IS_EMPTY(qt) || QT_LEAF_IS_FLAG(qt))
# Line 2150 | Line 2437 | mark_active_tris(argptr,root,qt,n)
2437    /* For each triangle in the set, set the which flag*/
2438    os = qtqueryset(qt);
2439  
2440 <  for (i = QT_SET_CNT(os), optr = QT_SET_PTR(os); i > 0; i--)
2440 >  for (i = QT_SET_CNT(os); i > 0; i--)
2441    {
2442 <    t_id = QT_SET_NEXT_ELEM(optr);
2443 <    /* Set the render flag */
2444 <     tri = SM_NTH_TRI(smMesh,t_id);
2445 <     if(!T_IS_VALID(tri))
2446 <            continue;
2447 <     SM_SET_NTH_T_ACTIVE(smMesh,t_id);
2448 <     /* Set the Active bits of the Vertices */
2449 <     S_SET_FLAG(T_NTH_V(tri,0));
2450 <     S_SET_FLAG(T_NTH_V(tri,1));
2451 <     S_SET_FLAG(T_NTH_V(tri,2));
2442 >    s_id = QT_SET_NEXT_ELEM(os);
2443 >    S_SET_FLAG(s_id);
2444 >
2445 >    tri_id = SM_NTH_VERT(smMesh,s_id);
2446 >    /* Set the active flag for all adjacent triangles */
2447 >     tri = SM_NTH_TRI(smMesh,tri_id);
2448 >     SM_SET_NTH_T_ACTIVE(smMesh,tri_id);
2449 >     while((t_id = smTri_next_ccw_nbr(smMesh,tri,s_id)) != tri_id)
2450 >     {
2451 >       tri = SM_NTH_TRI(smMesh,t_id);
2452 >       SM_SET_NTH_T_ACTIVE(smMesh,t_id);
2453 >     }
2454    }
2455 < }
2455 > }
2456  
2168
2457   smCull(sm,view,n)
2458   SM *sm;
2459   VIEW *view;
# Line 2187 | Line 2475 | smCull(sm,view,n)
2475          LRU counter: for use in discarding samples when out
2476          of space
2477          */
2478 <       F_FUNC(f) = mark_active_tris;
2478 >       F_FUNC(f) = mark_active_samples;
2479         smClear_flags(sm,T_ACTIVE_FLAG);
2480         /* Clear all of the active sample flags*/
2481         sClear_all_flags(SM_SAMP(sm));
# Line 2297 | Line 2585 | smCull(sm,view,n)
2585       }
2586  
2587   }
2588 +
2589 +
2590 +
2591 +
2592 +
2593 +
2594 +
2595 +
2596 +
2597 +
2598  
2599  
2600  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines