ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/sm.c
Revision: 3.3
Committed: Mon Aug 24 12:38:57 1998 UTC (25 years, 8 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.2: +66 -18 lines
Log Message:
optimized triangle addition/deletion during edge swapping

File Contents

# User Rev Content
1 gwlarson 3.1 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ SGI";
5     #endif
6    
7     /*
8     * sm.c
9     */
10     #include "standard.h"
11    
12     #include "object.h"
13    
14     #include "sm_list.h"
15     #include "sm_geom.h"
16     #include "sm.h"
17    
18    
19     SM *smMesh = NULL;
20     double smDist_sum=0;
21     int smNew_tri_cnt=0;
22    
23     #ifdef TEST_DRIVER
24     VIEW View = {0,{0,0,0},{0,0,-1},{0,1,0},60,60,0};
25     VIEW Current_View = {0,{0,0,0},{0,0,-1},{0,1,0},60,60,0};
26     int Pick_cnt =0;
27     int Pick_tri = -1,Picking = FALSE;
28     FVECT Pick_point[500],Pick_origin,Pick_dir;
29     FVECT Pick_v0[500], Pick_v1[500], Pick_v2[500];
30     FVECT P0,P1,P2;
31     FVECT FrustumNear[4],FrustumFar[4];
32     double dev_zmin=.01,dev_zmax=1000;
33     #endif
34    
35     smDir(sm,ps,id)
36     SM *sm;
37     FVECT ps;
38     int id;
39     {
40     FVECT p;
41    
42     VCOPY(p,SM_NTH_WV(sm,id));
43     point_on_sphere(ps,p,SM_VIEW_CENTER(sm));
44     }
45    
46     smClear_mesh(sm)
47     SM *sm;
48     {
49     /* Reset the triangle counters */
50     SM_TRI_CNT(sm) = 0;
51     SM_NUM_TRIS(sm) = 0;
52     SM_FREE_TRIS(sm) = -1;
53     }
54    
55     smClear_flags(sm,which)
56     SM *sm;
57     int which;
58     {
59     int i;
60    
61     if(which== -1)
62     for(i=0; i < T_FLAGS;i++)
63     bzero(SM_NTH_FLAGS(sm,i),T_TOTAL_FLAG_BYTES(SM_MAX_TRIS(sm)));
64     else
65     bzero(SM_NTH_FLAGS(sm,which),T_TOTAL_FLAG_BYTES(SM_MAX_TRIS(sm)));
66     }
67    
68     smClear_locator(sm)
69     SM *sm;
70     {
71     STREE *st;
72    
73     st = SM_LOCATOR(sm);
74    
75     stClear(st);
76     }
77    
78     smInit_locator(sm,center,base)
79     SM *sm;
80     FVECT center,base[4];
81     {
82     STREE *st;
83    
84     st = SM_LOCATOR(sm);
85    
86     stInit(st,center,base);
87    
88     }
89    
90     smClear(sm)
91     SM *sm;
92     {
93     smClear_samples(sm);
94     smClear_mesh(sm);
95     smClear_locator(sm);
96     }
97    
98     int
99     smAlloc_tris(sm,max_verts,max_tris)
100     SM *sm;
101     int max_verts,max_tris;
102     {
103     int i,nbytes,vbytes,fbytes;
104    
105     vbytes = max_verts*sizeof(VERT);
106     fbytes = T_TOTAL_FLAG_BYTES(max_tris);
107     nbytes = vbytes + max_tris*sizeof(TRI) +T_FLAGS*fbytes + 8;
108     for(i = 1024; nbytes > i; i <<= 1)
109     ;
110     /* check if casting works correctly */
111     max_tris = (i-vbytes-8)/(sizeof(TRI) + T_FLAG_BYTES);
112     fbytes = T_TOTAL_FLAG_BYTES(max_tris);
113    
114     SM_BASE(sm)=(char *)malloc(vbytes+max_tris*sizeof(TRI)+T_FLAGS*fbytes);
115    
116     if (SM_BASE(sm) == NULL)
117     return(0);
118    
119     SM_TRIS(sm) = (TRI *)SM_BASE(sm);
120     SM_VERTS(sm) = (VERT *)(SM_TRIS(sm) + max_tris);
121    
122     SM_NTH_FLAGS(sm,0) = (int4 *)(SM_VERTS(sm) + max_verts);
123     for(i=1; i < T_FLAGS; i++)
124     SM_NTH_FLAGS(sm,i) = (int4 *)(SM_NTH_FLAGS(sm,i-1)+fbytes/sizeof(int4));
125    
126     SM_MAX_VERTS(sm) = max_verts;
127     SM_MAX_TRIS(sm) = max_tris;
128    
129     smClear_mesh(sm);
130    
131     return(max_tris);
132     }
133    
134    
135    
136     int
137     smAlloc_locator(sm)
138     SM *sm;
139     {
140     STREE *st;
141    
142     st = SM_LOCATOR(sm);
143    
144     st = stAlloc(st);
145    
146     if(st)
147     return(TRUE);
148     else
149     return(FALSE);
150     }
151    
152     /* Initialize/clear global smL sample list for at least n samples */
153     smAlloc(max_samples)
154     register int max_samples;
155     {
156     unsigned nbytes;
157     register unsigned i;
158     int total_points;
159     int max_tris;
160    
161     /* If this is the first call, allocate sample,vertex and triangle lists */
162     if(!smMesh)
163     {
164     if(!(smMesh = (SM *)malloc(sizeof(SM))))
165     error(SYSTEM,"smAlloc():Unable to allocate memory");
166     bzero(smMesh,sizeof(SM));
167     }
168     else
169     { /* If existing structure: first deallocate */
170     if(SM_BASE(smMesh))
171     free(SM_BASE(smMesh));
172     if(SM_SAMP_BASE(smMesh))
173     free(SM_SAMP_BASE(smMesh));
174     }
175    
176     /* First allocate at least n samples + extra points:at least enough
177     necessary to form the BASE MESH- Default = 4;
178     */
179     max_samples = smAlloc_samples(smMesh, max_samples,SM_EXTRA_POINTS);
180    
181     total_points = max_samples + SM_EXTRA_POINTS;
182     max_tris = total_points*2;
183    
184     /* Now allocate space for mesh vertices and triangles */
185     max_tris = smAlloc_tris(smMesh, total_points, max_tris);
186    
187     /* Initialize the structure for point,triangle location.
188     */
189     smAlloc_locator(smMesh);
190    
191     }
192    
193    
194    
195     smInit_mesh(sm,vp)
196     SM *sm;
197     FVECT vp;
198     {
199    
200     /* NOTE: Should be elsewhere?*/
201     smDist_sum = 0;
202     smNew_tri_cnt = 0;
203    
204     VCOPY(SM_VIEW_CENTER(smMesh),vp);
205     smClear_locator(sm);
206     smInit_locator(sm,vp,0);
207     smClear_aux_samples(sm);
208     smClear_mesh(sm);
209     smCreate_base_mesh(sm,SM_DEFAULT);
210     }
211    
212     /*
213     * int
214     * smInit(n) : Initialize/clear data structures for n entries
215     * int n;
216     *
217     * This routine allocates/initializes the sample, mesh, and point-location
218     * structures for at least n samples.
219     * If n is <= 0, then clear data structures. Returns number samples
220     * actually allocated.
221     */
222    
223     int
224     smInit(n)
225     register int n;
226     {
227     int max_vertices;
228    
229     sleep(10);
230    
231     /* If n <=0, Just clear the existing structures */
232     if(n <= 0)
233     {
234     smClear(smMesh);
235     return(0);
236     }
237    
238     /* Total mesh vertices includes the sample points and the extra vertices
239     to form the base mesh
240     */
241     max_vertices = n + SM_EXTRA_POINTS;
242    
243     /* If the current mesh contains enough room, clear and return */
244     if(smMesh && max_vertices <= SM_MAX_VERTS(smMesh))
245     {
246     smClear(smMesh);
247     return(SM_MAX_SAMP(smMesh));
248     }
249     /* Otherwise- mesh must be allocated with the appropriate number of
250     samples
251     */
252     smAlloc(n);
253    
254     return(SM_MAX_SAMP(smMesh));
255     }
256    
257    
258     int
259     smLocator_apply_func(sm,v0,v1,v2,func,arg)
260     SM *sm;
261     FVECT v0,v1,v2;
262     int (*func)();
263     char *arg;
264     {
265     STREE *st;
266     char found;
267     FVECT p0,p1,p2;
268    
269     st = SM_LOCATOR(sm);
270    
271     point_on_sphere(p0,v0,SM_VIEW_CENTER(sm));
272     point_on_sphere(p1,v1,SM_VIEW_CENTER(sm));
273     point_on_sphere(p2,v2,SM_VIEW_CENTER(sm));
274    
275     found = stApply_to_tri_cells(st,p0,p1,p2,func,arg);
276    
277     return(found);
278     }
279    
280    
281     int
282     smLocator_add_tri(sm,t_id,v0_id,v1_id,v2_id)
283     SM *sm;
284     int t_id;
285     int v0_id,v1_id,v2_id;
286     {
287     STREE *st;
288     FVECT p0,p1,p2;
289    
290     st = SM_LOCATOR(sm);
291    
292     smDir(sm,p0,v0_id);
293     smDir(sm,p1,v1_id);
294     smDir(sm,p2,v2_id);
295    
296     stAdd_tri(st,t_id,p0,p1,p2);
297    
298     return(1);
299     }
300    
301     /* Add a triangle to the base array with vertices v1-v2-v3 */
302     int
303     smAdd_tri(sm, v0_id,v1_id,v2_id,tptr)
304     SM *sm;
305     int v0_id,v1_id,v2_id;
306     TRI **tptr;
307     {
308     int t_id;
309     TRI *t;
310    
311    
312     if(SM_TRI_CNT(sm)+1 > SM_MAX_TRIS(sm))
313     error(SYSTEM,"smAdd_tri():Too many triangles");
314    
315     /* Get the id for the next available triangle */
316     SM_FREE_TRI_ID(sm,t_id);
317     if(t_id == -1)
318     return(t_id);
319    
320     t = SM_NTH_TRI(sm,t_id);
321     T_CLEAR_NBRS(t);
322    
323     if(SM_BASE_ID(sm,v0_id) || SM_BASE_ID(sm,v1_id) || SM_BASE_ID(sm,v2_id))
324     {
325     smClear_tri_flags(sm,t_id);
326     SM_SET_NTH_T_BASE(sm,t_id);
327     }
328     else
329     {
330     SM_CLEAR_NTH_T_BASE(sm,t_id);
331     SM_SET_NTH_T_ACTIVE(sm,t_id);
332     SM_SET_NTH_T_LRU(sm,t_id);
333     SM_SET_NTH_T_NEW(sm,t_id);
334     SM_NUM_TRIS(sm)++;
335     smNew_tri_cnt++;
336     }
337     /* set the triangle vertex ids */
338     T_NTH_V(t,0) = v0_id;
339     T_NTH_V(t,1) = v1_id;
340     T_NTH_V(t,2) = v2_id;
341    
342     SM_NTH_VERT(sm,v0_id) = t_id;
343     SM_NTH_VERT(sm,v1_id) = t_id;
344     SM_NTH_VERT(sm,v2_id) = t_id;
345    
346     if(t)
347     *tptr = t;
348     /* return initialized triangle */
349     return(t_id);
350     }
351    
352     int
353     smClosest_vertex_in_tri(sm,v0_id,v1_id,v2_id,p,eps)
354     SM *sm;
355     int v0_id,v1_id,v2_id;
356     FVECT p;
357     double eps;
358     {
359     FVECT v;
360     double d,d0,d1,d2;
361     int closest = -1;
362    
363     if(v0_id != -1)
364     {
365     smDir(sm,v,v0_id);
366     d0 = DIST(p,v);
367     if(eps < 0 || d0 < eps)
368     {
369     closest = v0_id;
370     d = d0;
371     }
372     }
373     if(v1_id != -1)
374     {
375     smDir(sm,v,v1_id);
376     d1 = DIST(p,v);
377     if(closest== -1)
378     {
379     if(eps < 0 || d1 < eps)
380     {
381     closest = v1_id;
382     d = d1;
383     }
384     }
385     else
386     if(d1 < d2)
387     {
388     if((eps < 0) || d1 < eps)
389     {
390     closest = v1_id;
391     d = d1;
392     }
393     }
394     else
395     if((eps < 0) || d2 < eps)
396     {
397     closest = v2_id;
398     d = d2;
399     }
400     }
401     if(v2_id != -1)
402     {
403     smDir(sm,v,v2_id);
404     d2 = DIST(p,v);
405     if((eps < 0) || d2 < eps)
406     if(closest== -1 ||(d2 < d))
407     return(v2_id);
408     }
409     return(closest);
410     }
411    
412    
413     int
414     smClosest_vertex_in_w_tri(sm,v0_id,v1_id,v2_id,p)
415     SM *sm;
416     int v0_id,v1_id,v2_id;
417     FVECT p;
418     {
419     FVECT v;
420     double d,d0,d1,d2;
421     int closest;
422    
423     VCOPY(v,SM_NTH_WV(sm,v0_id));
424     d = d0 = DIST(p,v);
425     closest = v0_id;
426    
427     VCOPY(v,SM_NTH_WV(sm,v1_id));
428     d1 = DIST(p,v);
429     if(d1 < d2)
430     {
431     closest = v1_id;
432     d = d1;
433     }
434     VCOPY(v,SM_NTH_WV(sm,v2_id));
435     d2 = DIST(p,v);
436     if(d2 < d)
437     return(v2_id);
438     else
439     return(closest);
440     }
441    
442     void
443 gwlarson 3.3 smTris_swap_edge(sm,t_id,t1_id,e,e1,tn_id,tn1_id,add,del)
444 gwlarson 3.1 SM *sm;
445     int t_id,t1_id;
446     int e,e1;
447     int **tn_id,**tn1_id;
448 gwlarson 3.3 LIST **add,**del;
449 gwlarson 3.1
450     {
451     TRI *t,*t1;
452     TRI *ta,*tb;
453     int verts[3],enext,eprev,e1next,e1prev;
454     TRI *n;
455     FVECT p1,p2,p3;
456     int ta_id,tb_id;
457     /* swap diagonal (e relative to t, and e1 relative to t1)
458     defined by quadrilateral
459     formed by t,t1- swap for the opposite diagonal
460     */
461     t = SM_NTH_TRI(sm,t_id);
462     t1 = SM_NTH_TRI(sm,t1_id);
463     enext = (e+1)%3;
464     eprev = (e+2)%3;
465     e1next = (e1+1)%3;
466     e1prev = (e1+2)%3;
467     verts[e] = T_NTH_V(t,e);
468     verts[enext] = T_NTH_V(t1,e1prev);
469     verts[eprev] = T_NTH_V(t,eprev);
470     ta_id = smAdd_tri(sm,verts[0],verts[1],verts[2],&ta);
471 gwlarson 3.3 *add = push_data(*add,ta_id);
472 gwlarson 3.1
473     verts[e1] = T_NTH_V(t1,e1);
474     verts[e1next] = T_NTH_V(t,eprev);
475     verts[e1prev] = T_NTH_V(t1,e1prev);
476     tb_id = smAdd_tri(sm,verts[0],verts[1],verts[2],&tb);
477 gwlarson 3.3 *add = push_data(*add,tb_id);
478    
479 gwlarson 3.1 /* set the neighbors */
480     T_NTH_NBR(ta,e) = T_NTH_NBR(t1,e1next);
481     T_NTH_NBR(tb,e1) = T_NTH_NBR(t,enext);
482     T_NTH_NBR(ta,enext) = tb_id;
483     T_NTH_NBR(tb,e1next) = ta_id;
484     T_NTH_NBR(ta,eprev) = T_NTH_NBR(t,eprev);
485     T_NTH_NBR(tb,e1prev) = T_NTH_NBR(t1,e1prev);
486    
487     /* Reset neighbor pointers of original neighbors */
488     n = SM_NTH_TRI(sm,T_NTH_NBR(t,enext));
489     T_NTH_NBR(n,T_NTH_NBR_PTR(t_id,n)) = tb_id;
490     n = SM_NTH_TRI(sm,T_NTH_NBR(t,eprev));
491     T_NTH_NBR(n,T_NTH_NBR_PTR(t_id,n)) = ta_id;
492    
493     n = SM_NTH_TRI(sm,T_NTH_NBR(t1,e1next));
494     T_NTH_NBR(n,T_NTH_NBR_PTR(t1_id,n)) = ta_id;
495     n = SM_NTH_TRI(sm,T_NTH_NBR(t1,e1prev));
496     T_NTH_NBR(n,T_NTH_NBR_PTR(t1_id,n)) = tb_id;
497    
498     /* Delete two parent triangles */
499 gwlarson 3.3 *del = push_data(*del,t_id);
500     if(SM_IS_NTH_T_NEW(sm,t_id))
501     SM_CLEAR_NTH_T_NEW(sm,t_id);
502     else
503     SM_CLEAR_NTH_T_BASE(sm,t_id);
504     *del = push_data(*del,t1_id);
505     if(SM_IS_NTH_T_NEW(sm,t1_id))
506     SM_CLEAR_NTH_T_NEW(sm,t1_id);
507     else
508     SM_CLEAR_NTH_T_BASE(sm,t1_id);
509 gwlarson 3.1 *tn_id = ta_id;
510     *tn1_id = tb_id;
511     }
512    
513 gwlarson 3.3 smUpdate_locator(sm,add_list,del_list)
514     SM *sm;
515     LIST *add_list,*del_list;
516     {
517     int t_id;
518     TRI *t;
519     while(add_list)
520     {
521     t_id = pop_list(&add_list);
522     if(!SM_IS_NTH_T_NEW(sm,t_id) && !SM_IS_NTH_T_BASE(sm,t_id))
523     {
524     SM_SET_NTH_T_NEW(sm,t_id);
525     continue;
526     }
527     t = SM_NTH_TRI(sm,t_id);
528     smLocator_add_tri(sm,t_id,T_NTH_V(t,0),T_NTH_V(t,1),T_NTH_V(t,2));
529     }
530    
531     while(del_list)
532     {
533     t_id = pop_list(&del_list);
534     if(SM_IS_NTH_T_NEW(sm,t_id))
535     {
536     smDelete_tri(sm,t_id);
537     continue;
538     }
539     smLocator_remove_tri(sm,t_id);
540     smDelete_tri(sm,t_id);
541     }
542     }
543 gwlarson 3.1 /* MUST add check for constrained edges */
544 gwlarson 3.3 int
545 gwlarson 3.1 smFix_tris(sm,id,tlist)
546     SM *sm;
547     int id;
548     LIST *tlist;
549     {
550     TRI *t,*t_opp;
551 gwlarson 3.3 FVECT p,p1,p2,p3;
552     int e,e1,swapped = 0;
553 gwlarson 3.1 int t_id,t_opp_id;
554 gwlarson 3.3 LIST *add_list,*del_list;
555    
556    
557     add_list = del_list = NULL;
558     VSUB(p,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
559 gwlarson 3.1 while(tlist)
560     {
561 gwlarson 3.3 t_id = pop_list(&tlist);
562 gwlarson 3.1 t = SM_NTH_TRI(sm,t_id);
563     e = (T_WHICH_V(t,id)+1)%3;
564     t_opp_id = T_NTH_NBR(t,e);
565     t_opp = SM_NTH_TRI(sm,t_opp_id);
566    
567     smDir(sm,p1,T_NTH_V(t_opp,0));
568     smDir(sm,p2,T_NTH_V(t_opp,1));
569     smDir(sm,p3,T_NTH_V(t_opp,2));
570     if(point_in_cone(p,p1,p2,p3))
571     {
572     swapped = 1;
573     e1 = T_NTH_NBR_PTR(t_id,t_opp);
574     /* check list for t_opp and Remove if there */
575     remove_from_list(t_opp_id,&tlist);
576 gwlarson 3.3 smTris_swap_edge(sm,t_id,t_opp_id,e,e1,&t_id,&t_opp_id,
577     &add_list,&del_list);
578 gwlarson 3.1 tlist = push_data(tlist,t_id);
579     tlist = push_data(tlist,t_opp_id);
580     }
581     }
582 gwlarson 3.3 smUpdate_locator(sm,add_list,del_list);
583    
584 gwlarson 3.1 return(swapped);
585     }
586    
587     /* Give the vertex "id" and a triangle "t" that it belongs to- return the
588     next nbr in a counter clockwise order about vertex "id"
589     */
590     int
591     smTri_next_ccw_nbr(sm,t,id)
592     SM *sm;
593     TRI *t;
594     int id;
595     {
596     int t_id;
597     int tri;
598    
599     /* Want the edge for which "id" is the destination */
600     t_id = (T_WHICH_V(t,id)+ 2)% 3;
601     tri = T_NTH_NBR(t,t_id);
602     return(tri);
603     }
604    
605     void
606     smReplace_point(sm,tri,id,nid)
607     SM *sm;
608     TRI *tri;
609     int id,nid;
610     {
611     TRI *t;
612     int t_id;
613    
614     T_NTH_V(tri,T_WHICH_V(tri,id)) = nid;
615    
616     t_id = smTri_next_ccw_nbr(sm,t,nid);
617     while((t= SM_NTH_TRI(sm,t_id)) != tri)
618     {
619     T_NTH_V(t,T_WHICH_V(t,id)) = nid;
620     t_id = smTri_next_ccw_nbr(sm,t,nid);
621     }
622     }
623    
624    
625     smClear_tri_flags(sm,id)
626     SM *sm;
627     int id;
628     {
629     int i;
630    
631     for(i=0; i < T_FLAGS; i++)
632     SM_CLEAR_NTH_T_FLAG(sm,id,i);
633    
634     }
635    
636     /* Locate the point-id in the point location structure: */
637     int
638     smReplace_vertex(sm,c,dir,p,tri_id,snew_id,type,which)
639     SM *sm;
640     COLR c;
641     FVECT dir,p;
642     int tri_id,snew_id;
643     char type,which;
644     {
645     int n_id,s_id;
646     TRI *tri;
647    
648     tri = SM_NTH_TRI(sm,tri_id);
649     /* Get the sample that corresponds to the "which" vertex of "tri" */
650     /* NEED: have re-init that sets clock flag */
651     /* If this is a base-sample: create a new sample and replace
652     all references to the base sample with references to the
653     new sample
654     */
655     s_id = T_NTH_V(tri,which);
656     if(SM_BASE_ID(sm,s_id))
657     {
658     if(snew_id != -1)
659     n_id = smAdd_sample_point(sm,c,dir,p);
660     else
661     n_id = snew_id;
662     smReplace_point(sm,tri,s_id,n_id);
663     s_id = n_id;
664     }
665     else /* If the sample exists, reset the values */
666     /* NOTE: This test was based on the SPHERICAL coordinates
667     of the point: If we are doing a multiple-sample-per
668     SPHERICAL pixel: we will want to test for equality-
669     and do other processing: for now: SINGLE SAMPLE PER
670     PIXEL
671     */
672     /* NOTE: snew_id needs to be marked as invalid?*/
673     if(snew_id == -1)
674     smInit_sample(sm,s_id,c,dir,p);
675     else
676     smReset_sample(sm,s_id,snew_id);
677     return(s_id);
678     }
679    
680    
681     /* Locate the point-id in the point location structure: */
682     int
683     smInsert_point_in_tri(sm,c,dir,p,s_id,tri_id)
684     SM *sm;
685     COLR c;
686     FVECT dir,p;
687     int s_id,tri_id;
688     {
689     TRI *tri,*t0,*t1,*t2,*nbr;
690     int v0_id,v1_id,v2_id,n_id;
691     int t0_id,t1_id,t2_id;
692     LIST *tlist;
693     FVECT npt;
694    
695     if(s_id == SM_INVALID)
696     s_id = smAdd_sample_point(sm,c,dir,p);
697    
698     tri = SM_NTH_TRI(sm,tri_id);
699     v0_id = T_NTH_V(tri,0);
700     v1_id = T_NTH_V(tri,1);
701     v2_id = T_NTH_V(tri,2);
702    
703     n_id = -1;
704     if(SM_BASE_ID(sm,v0_id)||SM_BASE_ID(sm,v1_id)||SM_BASE_ID(sm,v2_id))
705     {
706     smDir(sm,npt,s_id);
707     /* Change to an add and delete */
708     t0_id = (SM_BASE_ID(sm,v0_id))?v0_id:-1;
709     t1_id = (SM_BASE_ID(sm,v1_id))?v1_id:-1;
710     t2_id = (SM_BASE_ID(sm,v2_id))?v2_id:-1;
711     n_id = smClosest_vertex_in_tri(sm,t0_id,t1_id,t2_id,npt,P_REPLACE_EPS);
712     }
713     t0_id = smAdd_tri(sm,s_id,v0_id,v1_id,&t0);
714 gwlarson 3.3 /* Add triangle to the locator */
715     smLocator_add_tri(sm,t0_id,s_id,v0_id,v1_id);
716    
717 gwlarson 3.1 t1_id = smAdd_tri(sm,s_id,v1_id,v2_id,&t1);
718 gwlarson 3.3 smLocator_add_tri(sm,t1_id,s_id,v1_id,v2_id);
719 gwlarson 3.1 t2_id = smAdd_tri(sm,s_id,v2_id,v0_id,&t2);
720 gwlarson 3.3 smLocator_add_tri(sm,t2_id,s_id,v2_id,v0_id);
721    
722 gwlarson 3.1 /* Set the neighbor pointers for the new tris */
723     T_NTH_NBR(t0,0) = t2_id;
724     T_NTH_NBR(t0,1) = T_NTH_NBR(tri,0);
725     T_NTH_NBR(t0,2) = t1_id;
726     T_NTH_NBR(t1,0) = t0_id;
727     T_NTH_NBR(t1,1) = T_NTH_NBR(tri,1);
728     T_NTH_NBR(t1,2) = t2_id;
729     T_NTH_NBR(t2,0) = t1_id;
730     T_NTH_NBR(t2,1) = T_NTH_NBR(tri,2);
731     T_NTH_NBR(t2,2) = t0_id;
732    
733     /* Reset the neigbor pointers for the neighbors of the original */
734     nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,0));
735     T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t0_id;
736     nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,1));
737     T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t1_id;
738     nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,2));
739     T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t2_id;
740    
741 gwlarson 3.3 smLocator_remove_tri(sm,tri_id);
742 gwlarson 3.1 smDelete_tri(sm,tri_id);
743    
744     /* Fix up the new triangles*/
745     tlist = push_data(NULL,t0_id);
746     tlist = push_data(tlist,t1_id);
747     tlist = push_data(tlist,t2_id);
748    
749     smFix_tris(sm,s_id,tlist);
750    
751     if(n_id != -1)
752     smDelete_point(sm,n_id);
753    
754     return(s_id);
755     }
756    
757    
758     int
759     smPointLocate(sm,pt,type,which,norm)
760     SM *sm;
761     FVECT pt;
762     char *type,*which;
763     char norm;
764     {
765     STREE *st;
766     int tri;
767     FVECT npt;
768    
769     st = SM_LOCATOR(sm);
770     if(norm)
771     {
772     point_on_sphere(npt,pt,SM_VIEW_CENTER(sm));
773     tri = stPoint_locate(st,npt,type,which);
774     }
775     else
776     tri = stPoint_locate(st,pt,type,which);
777     return(tri);
778     }
779    
780     QUADTREE
781 gwlarson 3.2 smPointLocateCell(sm,pt,type,which,norm)
782 gwlarson 3.1 SM *sm;
783 gwlarson 3.2 FVECT pt;
784 gwlarson 3.1 char *type,*which;
785     char norm;
786     {
787     STREE *st;
788     QUADTREE qt;
789     FVECT npt;
790    
791     st = SM_LOCATOR(sm);
792     if(norm)
793     {
794     point_on_sphere(npt,pt,SM_VIEW_CENTER(sm));
795    
796 gwlarson 3.2 qt = stPoint_locate_cell(st,npt,type,which);
797 gwlarson 3.1 }
798     else
799 gwlarson 3.2 qt = stPoint_locate_cell(st,pt,type,which);
800 gwlarson 3.1
801     return(qt);
802     }
803    
804     int
805     smAdd_sample_to_mesh(sm,c,dir,pt,s_id)
806     SM *sm;
807     COLR c;
808     FVECT dir,pt;
809     int s_id;
810     {
811     int t_id;
812     char type,which;
813     double d;
814     FVECT p;
815    
816     /* If new, foreground pt */
817     if(pt)
818     {
819     /* NOTE: This should be elsewhere! */
820     d = DIST(pt,SM_VIEW_CENTER(smMesh));
821     smDist_sum += 1.0/d;
822     /************************************/
823     t_id = smPointLocate(smMesh,pt,&type,&which,TRUE);
824     if(type==GT_FACE)
825     s_id = smInsert_point_in_tri(smMesh,c,dir,pt,s_id,t_id);
826     else
827     if(type==GT_VERTEX)
828     s_id = smReplace_vertex(smMesh,c,dir,pt,t_id,s_id,type,which);
829     #ifdef DEBUG
830     else
831     eputs("smAdd_sample_to_mesh(): unrecognized type\n");
832     #endif
833     }
834     else if(s_id != -1)
835     {
836     VCOPY(p,SM_NTH_WV(sm,s_id));
837     if(SM_NTH_W_DIR(sm,s_id) != -1)
838     {
839     /* NOTE: This should be elsewhere! */
840     d = DIST(p,SM_VIEW_CENTER(smMesh));
841     smDist_sum += 1.0/d;
842     /************************************/
843     }
844     t_id = smPointLocate(smMesh,p,&type,&which,TRUE);
845     if(type==GT_FACE)
846     s_id = smInsert_point_in_tri(smMesh,c,dir,p,s_id,t_id);
847     else
848     if(type==GT_VERTEX)
849     s_id = smReplace_vertex(smMesh,c,dir,p,t_id,s_id,type,which);
850     #ifdef DEBUG
851     else
852     eputs("smAdd_sample_to_mesh(): unrecognized type\n");
853     #endif
854     }
855     /* Is a BG(sky point) */
856     else
857     {
858     t_id = smPointLocate(smMesh,dir,&type,&which,FALSE);
859     if(type==GT_FACE)
860     s_id = smInsert_point_in_tri(smMesh,c,dir,NULL,s_id,t_id);
861     else
862     if(type==GT_VERTEX)
863     s_id = smReplace_vertex(smMesh,c,dir,NULL,t_id,s_id,type,which);
864     #ifdef DEBUG
865     else
866     eputs("smAdd_sample_to_mesh(): unrecognized type\n");
867     #endif
868     }
869     return(s_id);
870     }
871    
872     /*
873     * int
874     * smNewSamp(c, dir, p) : register new sample point and return index
875     * COLR c; : pixel color (RGBE)
876     * FVECT dir; : ray direction vector
877     * FVECT p; : world intersection point
878     *
879     * Add new sample point to data structures, removing old values as necessary.
880     * New sample representation will be output in next call to smUpdate().
881     * If the point is a sky point: then v=NULL
882     */
883     int
884     smNewSamp(c,dir,p)
885     COLR c;
886     FVECT dir;
887     FVECT p;
888    
889     {
890     int s_id;
891    
892     /* First check if this the first sample: if so initialize mesh */
893     if(SM_NUM_SAMP(smMesh) == 0)
894     #ifdef TEST_DRIVER
895     smInit_mesh(smMesh,View.vp);
896     #else
897     smInit_mesh(smMesh,odev.v.vp);
898     #endif
899     s_id = smAdd_sample_to_mesh(smMesh,c,dir,p,-1);
900     #if 0
901     {
902     char buff[100];
903     sprintf(buff,"Added sample %d\n",s_id);
904     eputs(buff);
905     }
906     #endif
907     return(s_id);
908    
909     }
910     /*
911     * int
912     * smFindsamp(orig, dir): intersect ray with 3D rep. and find closest sample
913     * FVECT orig, dir;
914     *
915     * Find the closest sample to the given ray. Return -1 on failure.
916     */
917    
918     /*
919     * smClean() : display has been wiped clean
920     *
921     * Called after display has been effectively cleared, meaning that all
922     * geometry must be resent down the pipeline in the next call to smUpdate().
923     */
924    
925    
926     /*
927     * smUpdate(vp, qua) : update OpenGL output geometry for view vp
928     * VIEW *vp; : desired view
929     * int qua; : quality level (percentage on linear time scale)
930     *
931     * Draw new geometric representation using OpenGL calls. Assume that the
932     * view has already been set up and the correct frame buffer has been
933     * selected for drawing. The quality level is on a linear scale, where 100%
934     * is full (final) quality. It is not necessary to redraw geometry that has
935     * been output since the last call to smClean().
936     */
937    
938    
939     int
940     smClear_vert(sm,id)
941     SM *sm;
942     int id;
943     {
944     if(SM_INVALID_POINT_ID(sm,id))
945     return(FALSE);
946    
947     SM_NTH_VERT(sm,id) = SM_INVALID;
948    
949     return(TRUE);
950     }
951    
952     int
953     smAdd_base_vertex(sm,v,d)
954     SM *sm;
955     FVECT v,d;
956     {
957     int id;
958    
959     /* First add coordinate to the sample array */
960     id = smAdd_aux_point(sm,v,d);
961     if(id == -1)
962     return(SM_INVALID);
963     /* Initialize triangle pointer to -1 */
964     smClear_vert(sm,id);
965     return(id);
966     }
967    
968    
969    
970     /* Initialize a the point location DAG based on a 6 triangle tesselation
971     of the unit sphere centered on the view center. The DAG structure
972     contains 6 roots: one for each initial base triangle
973     */
974    
975     int
976     smCreate_base_mesh(sm,type)
977     SM *sm;
978     int type;
979     {
980     int i,id;
981     int p[4],ids[4];
982     int v0_id,v1_id,v2_id;
983     TRI *tris[4];
984     FVECT d,pt,cntr;
985    
986     /* First insert the base vertices into the sample point array */
987    
988     for(i=0; i < 4; i++)
989     {
990     VADD(cntr,stDefault_base[i],SM_VIEW_CENTER(sm));
991     point_on_sphere(d,cntr,SM_VIEW_CENTER(sm));
992     id = smAdd_base_vertex(sm,cntr,d);
993     /* test to make sure vertex allocated */
994     if(id != -1)
995     p[i] = id;
996     else
997     return(0);
998     }
999     /* Create the base triangles */
1000     for(i=0;i < 4; i++)
1001     {
1002     v0_id = p[stTri_verts[i][0]];
1003     v1_id = p[stTri_verts[i][1]];
1004     v2_id = p[stTri_verts[i][2]];
1005     if((ids[i] = smAdd_tri(sm, v0_id,v1_id,v2_id,&(tris[i])))== -1)
1006     return(0);
1007 gwlarson 3.3 smLocator_add_tri(sm,ids[i],v0_id,v1_id,v2_id);
1008 gwlarson 3.1 }
1009     /* Set neighbors */
1010    
1011     T_NTH_NBR(tris[0],0) = ids[3];
1012     T_NTH_NBR(tris[0],1) = ids[2];
1013     T_NTH_NBR(tris[0],2) = ids[1];
1014    
1015     T_NTH_NBR(tris[1],0) = ids[3];
1016     T_NTH_NBR(tris[1],1) = ids[0];
1017     T_NTH_NBR(tris[1],2) = ids[2];
1018    
1019     T_NTH_NBR(tris[2],0) = ids[3];
1020     T_NTH_NBR(tris[2],1) = ids[1];
1021     T_NTH_NBR(tris[2],2) = ids[0];
1022    
1023     T_NTH_NBR(tris[3],0) = ids[1];
1024     T_NTH_NBR(tris[3],1) = ids[2];
1025     T_NTH_NBR(tris[3],2) = ids[0];
1026     return(1);
1027    
1028     }
1029    
1030    
1031     int
1032     smNext_tri_flag_set(sm,i,which,b)
1033     SM *sm;
1034     int i,which;
1035     char b;
1036     {
1037    
1038     for(; i < SM_TRI_CNT(sm);i++)
1039     {
1040     if(!SM_IS_NTH_T_FLAG(sm,i,which))
1041     continue;
1042    
1043     if(!b)
1044     break;
1045     if((b==1) && !SM_BG_TRI(sm,i))
1046     break;
1047     if((b==2) && SM_BG_TRI(sm,i))
1048     break;
1049     }
1050    
1051     return(i);
1052     }
1053    
1054    
1055     int
1056     smNext_valid_tri(sm,i)
1057     SM *sm;
1058     int i;
1059     {
1060    
1061     while( i < SM_TRI_CNT(sm) && !T_IS_VALID(SM_NTH_TRI(sm,i)))
1062     i++;
1063    
1064     return(i);
1065     }
1066    
1067    
1068    
1069     qtTri_verts_from_id(t_id,v0,v1,v2)
1070     int t_id;
1071     FVECT v0,v1,v2;
1072     {
1073     TRI *t;
1074     int v0_id,v1_id,v2_id;
1075    
1076     t = SM_NTH_TRI(smMesh,t_id);
1077     v0_id = T_NTH_V(t,0);
1078     v1_id = T_NTH_V(t,1);
1079     v2_id = T_NTH_V(t,2);
1080    
1081     smDir(smMesh,v0,v0_id);
1082     smDir(smMesh,v1,v1_id);
1083     smDir(smMesh,v2,v2_id);
1084    
1085     }
1086    
1087    
1088     int
1089     smIntersectTriSet(sm,t_set,orig,dir,pt)
1090     SM *sm;
1091     OBJECT *t_set;
1092     FVECT orig,dir,pt;
1093     {
1094     OBJECT *optr;
1095     int i,t_id,v_id;
1096     TRI *tri;
1097     FVECT p0,p1,p2;
1098     char type,which;
1099     int p0_id,p1_id,p2_id;
1100    
1101     for(optr = QT_SET_PTR(t_set),i = QT_SET_CNT(t_set); i > 0; i--)
1102     {
1103     t_id = QT_SET_NEXT_ELEM(optr);
1104     tri = SM_NTH_TRI(sm,t_id);
1105     p0_id = T_NTH_V(tri,0);
1106     p1_id = T_NTH_V(tri,1);
1107     p2_id = T_NTH_V(tri,2);
1108     VCOPY(p0,SM_NTH_WV(sm,p0_id));
1109     VCOPY(p1,SM_NTH_WV(sm,p1_id));
1110     VCOPY(p2,SM_NTH_WV(sm,p2_id));
1111     if(type = ray_intersect_tri(orig,dir,p0,p1,p2,pt,&which))
1112     {
1113     if(type==GT_VERTEX)
1114     return(T_NTH_V(tri,which));
1115     v_id = smClosest_vertex_in_w_tri(sm,p0_id,p1_id,p2_id,pt);
1116     return(v_id);
1117     }
1118     }
1119     return(-1);
1120     }
1121    
1122     /*
1123     * int
1124     * smTraceRay(SM *sm,FVECT orig, FVECT dir,FVECT v0,FVECT v1,FVECT v2,FVECT r)
1125     *
1126     * Intersect the ray with triangle v0v1v2, return intersection point in r
1127     *
1128     * Assumes orig,v0,v1,v2 are in spherical coordinates, and orig is
1129     * unit
1130     */
1131     char
1132     smTraceRay(sm,orig,dir,v0,v1,v2,r)
1133     SM *sm;
1134     FVECT orig,dir;
1135     FVECT v0,v1,v2;
1136     FVECT r;
1137     {
1138     FVECT n,p[3],d;
1139     double pt[3],r_eps;
1140     char i;
1141     int which;
1142    
1143     /* Find the plane equation for the triangle defined by the edge v0v1 and
1144     the view center*/
1145     VCROSS(n,v0,v1);
1146     /* Intersect the ray with this plane */
1147     i = intersect_ray_plane(orig,dir,n,0.0,&(pt[0]),p[0]);
1148     if(i)
1149     which = 0;
1150     else
1151     which = -1;
1152    
1153     VCROSS(n,v1,v2);
1154     i = intersect_ray_plane(orig,dir,n,0.0,&(pt[1]),p[1]);
1155     if(i)
1156     if((which==-1) || pt[1] < pt[0])
1157     which = 1;
1158    
1159     VCROSS(n,v2,v0);
1160     i = intersect_ray_plane(orig,dir,n,0.0,&(pt[2]),p[2]);
1161     if(i)
1162     if((which==-1) || pt[2] < pt[which])
1163     which = 2;
1164    
1165     if(which != -1)
1166     {
1167     /* Return point that is K*eps along projection of the ray on
1168     the sphere to push intersection point p[which] into next cell
1169     */
1170     normalize(p[which]);
1171     /* Calculate the ray perpendicular to the intersection point: approx
1172     the direction moved along the sphere a distance of K*epsilon*/
1173     r_eps = -DOT(p[which],dir);
1174     VSUM(n,dir,p[which],r_eps);
1175     /* Calculate the length along ray p[which]-dir needed to move to
1176     cause a move along the sphere of k*epsilon
1177     */
1178     r_eps = DOT(n,dir);
1179     VSUM(r,p[which],dir,(20*FTINY)/r_eps);
1180     normalize(r);
1181     return(TRUE);
1182     }
1183     else
1184     {
1185     /* Unable to find intersection: move along ray and try again */
1186     r_eps = -DOT(orig,dir);
1187     VSUM(n,dir,orig,r_eps);
1188     r_eps = DOT(n,dir);
1189     VSUM(r,orig,dir,(20*FTINY)/r_eps);
1190     normalize(r);
1191     #ifdef DEBUG
1192     eputs("smTraceRay:Ray does not intersect triangle");
1193     #endif
1194     return(FALSE);
1195     }
1196     }
1197    
1198    
1199     /*
1200     * int
1201     * smFindSamp(FVECT orig, FVECT dir)
1202     *
1203     * Find the closest sample to the given ray. Returns sample id, -1 on failure.
1204     * "dir" is assumed to be normalized
1205     */
1206     int
1207     smFindSamp(orig,dir)
1208     FVECT orig,dir;
1209     {
1210     FVECT r,v0,v1,v2,a,b,p;
1211     OBJECT os[MAXCSET+1],t_set[MAXSET+1];
1212     QUADTREE qt;
1213     int s_id;
1214     double d;
1215    
1216     /* r is the normalized vector from the view center to the current
1217     * ray point ( starting with "orig"). Find the cell that r falls in,
1218     * and test the ray against all triangles stored in the cell. If
1219     * the test fails, trace the projection of the ray across to the
1220     * next cell it intersects: iterate until either an intersection
1221     * is found, or the projection ray is // to the direction. The sample
1222     * corresponding to the triangle vertex closest to the intersection
1223     * point is returned.
1224     */
1225    
1226     /* First test if "orig" coincides with the View_center or if "dir" is
1227     parallel to r formed by projecting "orig" on the sphere. In
1228     either case, do a single test against the cell containing the
1229     intersection of "dir" and the sphere
1230     */
1231     point_on_sphere(b,orig,SM_VIEW_CENTER(smMesh));
1232     d = -DOT(b,dir);
1233     if(EQUAL_VEC3(orig,SM_VIEW_CENTER(smMesh)) || EQUAL(fabs(d),1.0))
1234     {
1235 gwlarson 3.2 qt = smPointLocateCell(smMesh,dir,NULL,NULL,FALSE);
1236 gwlarson 3.1 /* Test triangles in the set for intersection with Ray:returns
1237     first found
1238     */
1239     qtgetset(t_set,qt);
1240     s_id = smIntersectTriSet(smMesh,t_set,orig,dir,p);
1241     #ifdef TEST_DRIVER
1242     VCOPY(Pick_point[0],p);
1243     #endif
1244     return(s_id);
1245     }
1246     else
1247     {
1248     /* Starting with orig, Walk along projection of ray onto sphere */
1249     point_on_sphere(r,orig,SM_VIEW_CENTER(smMesh));
1250 gwlarson 3.2 qt = smPointLocateCell(smMesh,r,NULL,NULL,FALSE);
1251 gwlarson 3.1 qtgetset(t_set,qt);
1252     /* os will contain all triangles seen thus far */
1253     setcopy(os,t_set);
1254    
1255     /* Calculate ray perpendicular to dir: when projection ray is // to dir,
1256     the dot product will become negative.
1257     */
1258     VSUM(a,b,dir,d);
1259     d = DOT(a,b);
1260     while(d > 0)
1261     {
1262     s_id = smIntersectTriSet(smMesh,t_set,orig,dir,p);
1263     #ifdef TEST_DRIVER
1264     VCOPY(Pick_point[0],p);
1265     #endif
1266     if(s_id != EMPTY)
1267     return(s_id);
1268     /* Find next cell that projection of ray intersects */
1269     smTraceRay(smMesh,r,dir,v0,v1,v2,r);
1270 gwlarson 3.2 qt = smPointLocateCell(smMesh,r,NULL,NULL,FALSE);
1271 gwlarson 3.1 qtgetset(t_set,qt);
1272     /* Check triangles in set against those seen so far(os):only
1273     check new triangles for intersection (t_set')
1274     */
1275     check_set(t_set,os);
1276     d = DOT(a,r);
1277     }
1278     }
1279     #ifdef DEBUG
1280     eputs("smFindSamp():Pick Ray did not intersect mesh");
1281     #endif
1282     return(EMPTY);
1283     }
1284    
1285    
1286     smRebuild_mesh(sm,vptr)
1287     SM *sm;
1288     VIEW *vptr;
1289     {
1290     int i;
1291     FVECT dir;
1292     COLR c;
1293     FVECT p,ov;
1294    
1295     /* Clear the mesh- and rebuild using the current sample array */
1296     #ifdef TEST_DRIVER
1297     View = *vptr;
1298     #endif
1299    
1300     VSUB(ov,vptr->vp,SM_VIEW_CENTER(sm));
1301     smInit_mesh(sm,vptr->vp);
1302    
1303     SM_FOR_ALL_SAMPLES(sm,i)
1304     {
1305     if(SM_NTH_W_DIR(sm,i)==-1)
1306     VADD(SM_NTH_WV(sm,i),SM_NTH_WV(sm,i),ov);
1307     smAdd_sample_to_mesh(sm,NULL,NULL,NULL,i);
1308     }
1309     }
1310    
1311