ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/sm.c
Revision: 3.8
Committed: Tue Oct 6 18:16:54 1998 UTC (25 years, 6 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.7: +841 -718 lines
Log Message:
new triangulate routine
added smTestSample to check for occlusion
added frustum culling before rebuild
changed base quadtree to use octahedron and created new point locate
added "sample active" flags and implemented LRU replacement
started handling case of too many triangles
set sizes are now unbounded
changed all quadtree pointers to quadtrees

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 gwlarson 3.8 #include "sm_flag.h"
12 gwlarson 3.1 #include "sm_list.h"
13     #include "sm_geom.h"
14 gwlarson 3.8 #include "sm_qtree.h"
15     #include "sm_stree.h"
16 gwlarson 3.1 #include "sm.h"
17    
18    
19     SM *smMesh = NULL;
20     double smDist_sum=0;
21     int smNew_tri_cnt=0;
22 gwlarson 3.8 double smMinSampDiff = 1e-4;
23 gwlarson 3.1
24 gwlarson 3.8 static FVECT smDefault_base[4] = { {SQRT3_INV, SQRT3_INV, SQRT3_INV},
25     {-SQRT3_INV, -SQRT3_INV, SQRT3_INV},
26     {-SQRT3_INV, SQRT3_INV, -SQRT3_INV},
27     {SQRT3_INV, -SQRT3_INV, -SQRT3_INV}};
28     static int smTri_verts[4][3] = { {2,1,0},{3,2,0},{1,3,0},{2,3,1}};
29    
30 gwlarson 3.5 static int smBase_nbrs[4][3] = { {3,2,1},{3,0,2},{3,1,0},{1,2,0}};
31    
32 gwlarson 3.1 #ifdef TEST_DRIVER
33     VIEW Current_View = {0,{0,0,0},{0,0,-1},{0,1,0},60,60,0};
34     int Pick_cnt =0;
35 gwlarson 3.5 int Pick_tri = -1,Picking = FALSE,Pick_samp=-1;
36 gwlarson 3.1 FVECT Pick_point[500],Pick_origin,Pick_dir;
37     FVECT Pick_v0[500], Pick_v1[500], Pick_v2[500];
38     FVECT P0,P1,P2;
39     FVECT FrustumNear[4],FrustumFar[4];
40     double dev_zmin=.01,dev_zmax=1000;
41     #endif
42    
43     smDir(sm,ps,id)
44     SM *sm;
45     FVECT ps;
46     int id;
47     {
48 gwlarson 3.8 VSUB(ps,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
49     normalize(ps);
50 gwlarson 3.1 }
51 gwlarson 3.8
52 gwlarson 3.7 smDir_in_cone(sm,ps,id)
53     SM *sm;
54     FVECT ps;
55     int id;
56     {
57    
58 gwlarson 3.8 VSUB(ps,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
59     normalize(ps);
60 gwlarson 3.7 }
61 gwlarson 3.1
62     smClear_flags(sm,which)
63     SM *sm;
64     int which;
65     {
66     int i;
67    
68     if(which== -1)
69     for(i=0; i < T_FLAGS;i++)
70 gwlarson 3.8 bzero(SM_NTH_FLAGS(sm,i),FLAG_BYTES(SM_MAX_TRIS(sm)));
71 gwlarson 3.1 else
72 gwlarson 3.8 bzero(SM_NTH_FLAGS(sm,which),FLAG_BYTES(SM_MAX_TRIS(sm)));
73 gwlarson 3.1 }
74    
75 gwlarson 3.8 /* Given an allocated mesh- initialize */
76     smInit_mesh(sm)
77     SM *sm;
78 gwlarson 3.1 {
79 gwlarson 3.8 /* Reset the triangle counters */
80     SM_NUM_TRI(sm) = 0;
81     SM_SAMPLE_TRIS(sm) = 0;
82     SM_FREE_TRIS(sm) = -1;
83     smClear_flags(sm,-1);
84 gwlarson 3.1 }
85    
86    
87 gwlarson 3.8 /* Clear the SM for reuse: free any extra allocated memory:does initialization
88     as well
89     */
90 gwlarson 3.1 smClear(sm)
91     SM *sm;
92     {
93     smClear_samples(sm);
94     smClear_locator(sm);
95 gwlarson 3.8 smInit_mesh(sm);
96 gwlarson 3.1 }
97    
98 gwlarson 3.8 static int realloc_cnt=0;
99    
100 gwlarson 3.1 int
101 gwlarson 3.8 smRealloc_mesh(sm)
102 gwlarson 3.1 SM *sm;
103 gwlarson 3.8 {
104     int fbytes,i,max_tris,max_verts;
105    
106     if(realloc_cnt <=0)
107     realloc_cnt = SM_MAX_TRIS(sm);
108    
109     max_tris = SM_MAX_TRIS(sm) + realloc_cnt;
110     fbytes = FLAG_BYTES(max_tris);
111    
112     if(!(SM_TRIS(sm) = (TRI *)realloc(SM_TRIS(sm),max_tris*sizeof(TRI))))
113     goto memerr;
114    
115     for(i=0; i< T_FLAGS; i++)
116     if(!(SM_NTH_FLAGS(sm,i)=(int4 *)realloc((char *)SM_NTH_FLAGS(sm,i),fbytes)))
117     goto memerr;
118    
119     SM_MAX_TRIS(sm) = max_tris;
120     return(max_tris);
121    
122     memerr:
123     error(SYSTEM, "out of memory in smRealloc_mesh()");
124     }
125    
126     /* Allocate and initialize a new mesh with max_verts and max_tris */
127     int
128     smAlloc_mesh(sm,max_verts,max_tris)
129     SM *sm;
130 gwlarson 3.1 int max_verts,max_tris;
131     {
132 gwlarson 3.8 int fbytes,i;
133 gwlarson 3.1
134 gwlarson 3.8 fbytes = FLAG_BYTES(max_tris);
135 gwlarson 3.1
136 gwlarson 3.8 if(!(SM_TRIS(sm) = (TRI *)realloc(NULL,max_tris*sizeof(TRI))))
137     goto memerr;
138 gwlarson 3.1
139 gwlarson 3.8 if(!(SM_VERTS(sm) = (VERT *)realloc(NULL,max_verts*sizeof(VERT))))
140     goto memerr;
141 gwlarson 3.1
142 gwlarson 3.8 for(i=0; i< T_FLAGS; i++)
143     if(!(SM_NTH_FLAGS(sm,i)=(int4 *)realloc(NULL,fbytes)))
144     goto memerr;
145    
146 gwlarson 3.1 SM_MAX_VERTS(sm) = max_verts;
147     SM_MAX_TRIS(sm) = max_tris;
148    
149 gwlarson 3.8 realloc_cnt = max_verts >> 1;
150 gwlarson 3.1
151 gwlarson 3.8 smInit_mesh(sm);
152    
153 gwlarson 3.1 return(max_tris);
154 gwlarson 3.8 memerr:
155     error(SYSTEM, "out of memory in smAlloc_mesh()");
156 gwlarson 3.1 }
157    
158    
159     int
160 gwlarson 3.8 smAlloc_tri(sm)
161 gwlarson 3.1 SM *sm;
162     {
163 gwlarson 3.8 int id;
164    
165     /* First check if there are any tris on the free list */
166     /* Otherwise, have we reached the end of the list yet?*/
167     if(SM_NUM_TRI(sm) < SM_MAX_TRIS(sm))
168     return(SM_NUM_TRI(sm)++);
169    
170     if((id = SM_FREE_TRIS(sm)) != -1)
171     {
172     SM_FREE_TRIS(sm) = T_NEXT_FREE(SM_NTH_TRI(sm,id));
173     return(id);
174     }
175    
176     /*Else: realloc */
177     smRealloc_mesh(sm);
178     return(SM_NUM_TRI(sm)++);
179 gwlarson 3.1 }
180    
181 gwlarson 3.8 smFree_mesh(sm)
182     SM *sm;
183     {
184     int i;
185    
186     if(SM_TRIS(sm))
187     free(SM_TRIS(sm));
188     if(SM_VERTS(sm))
189     free(SM_VERTS(sm));
190     for(i=0; i< T_FLAGS; i++)
191     if(SM_NTH_FLAGS(sm,i))
192     free(SM_NTH_FLAGS(sm,i));
193     }
194    
195    
196 gwlarson 3.1 /* Initialize/clear global smL sample list for at least n samples */
197     smAlloc(max_samples)
198     register int max_samples;
199     {
200     unsigned nbytes;
201     register unsigned i;
202     int total_points;
203 gwlarson 3.8 int max_tris,n;
204 gwlarson 3.1
205 gwlarson 3.8 n = max_samples;
206 gwlarson 3.1 /* If this is the first call, allocate sample,vertex and triangle lists */
207     if(!smMesh)
208     {
209     if(!(smMesh = (SM *)malloc(sizeof(SM))))
210 gwlarson 3.8 error(SYSTEM,"smAlloc():Unable to allocate memory\n");
211 gwlarson 3.1 bzero(smMesh,sizeof(SM));
212     }
213     else
214     { /* If existing structure: first deallocate */
215 gwlarson 3.8 smFree_mesh(smMesh);
216     smFree_samples(smMesh);
217     smFree_locator(smMesh);
218     }
219 gwlarson 3.1
220     /* First allocate at least n samples + extra points:at least enough
221     necessary to form the BASE MESH- Default = 4;
222     */
223 gwlarson 3.8 SM_SAMP(smMesh) = sAlloc(&n,SM_EXTRA_POINTS);
224 gwlarson 3.1
225 gwlarson 3.8 total_points = n + SM_EXTRA_POINTS;
226 gwlarson 3.1
227 gwlarson 3.8 max_tris = total_points*4;
228 gwlarson 3.1 /* Now allocate space for mesh vertices and triangles */
229 gwlarson 3.8 max_tris = smAlloc_mesh(smMesh, total_points, max_tris);
230 gwlarson 3.1
231     /* Initialize the structure for point,triangle location.
232     */
233     smAlloc_locator(smMesh);
234     }
235    
236    
237    
238 gwlarson 3.8 smInit_sm(sm,vp)
239 gwlarson 3.1 SM *sm;
240     FVECT vp;
241     {
242    
243     smDist_sum = 0;
244     smNew_tri_cnt = 0;
245    
246 gwlarson 3.8 VCOPY(SM_VIEW_CENTER(sm),vp);
247     smInit_locator(sm,vp);
248     smInit_samples(sm);
249     smInit_mesh(sm);
250 gwlarson 3.1 smCreate_base_mesh(sm,SM_DEFAULT);
251     }
252    
253     /*
254     * int
255     * smInit(n) : Initialize/clear data structures for n entries
256     * int n;
257     *
258     * This routine allocates/initializes the sample, mesh, and point-location
259     * structures for at least n samples.
260     * If n is <= 0, then clear data structures. Returns number samples
261     * actually allocated.
262     */
263    
264     int
265     smInit(n)
266     register int n;
267     {
268     int max_vertices;
269    
270 gwlarson 3.8
271 gwlarson 3.1 /* If n <=0, Just clear the existing structures */
272     if(n <= 0)
273     {
274     smClear(smMesh);
275     return(0);
276     }
277    
278     /* Total mesh vertices includes the sample points and the extra vertices
279     to form the base mesh
280     */
281     max_vertices = n + SM_EXTRA_POINTS;
282    
283     /* If the current mesh contains enough room, clear and return */
284 gwlarson 3.8 if(smMesh && (max_vertices <= SM_MAX_VERTS(smMesh)) && SM_MAX_SAMP(smMesh) <=
285     n && SM_MAX_POINTS(smMesh) <= max_vertices)
286 gwlarson 3.1 {
287     smClear(smMesh);
288     return(SM_MAX_SAMP(smMesh));
289     }
290     /* Otherwise- mesh must be allocated with the appropriate number of
291     samples
292     */
293     smAlloc(n);
294    
295     return(SM_MAX_SAMP(smMesh));
296     }
297    
298    
299 gwlarson 3.8 smLocator_apply_func(sm,v0,v1,v2,edge_func,tri_func,argptr)
300 gwlarson 3.7 SM *sm;
301     FVECT v0,v1,v2;
302     int (*edge_func)();
303 gwlarson 3.8 int (*tri_func)();
304     int *argptr;
305 gwlarson 3.1 {
306     STREE *st;
307     FVECT p0,p1,p2;
308    
309     st = SM_LOCATOR(sm);
310    
311 gwlarson 3.5 VSUB(p0,v0,SM_VIEW_CENTER(sm));
312     VSUB(p1,v1,SM_VIEW_CENTER(sm));
313     VSUB(p2,v2,SM_VIEW_CENTER(sm));
314 gwlarson 3.1
315 gwlarson 3.8 stApply_to_tri(st,p0,p1,p2,edge_func,tri_func,argptr);
316 gwlarson 3.1
317     }
318    
319 gwlarson 3.8 QUADTREE
320     delete_tri_set(qt,optr,del_set)
321     QUADTREE qt;
322     OBJECT *optr,*del_set;
323     {
324 gwlarson 3.1
325 gwlarson 3.8 int set_size,i,t_id;
326     OBJECT *rptr,r_set[QT_MAXSET+1];
327     OBJECT *tptr,t_set[QT_MAXSET+1],*sptr;
328    
329     /* First need to check if set size larger than QT_MAXSET: if so
330     need to allocate larger array
331     */
332     if((set_size = MAX(QT_SET_CNT(optr),QT_SET_CNT(del_set))) >QT_MAXSET)
333     rptr = (OBJECT *)malloc((set_size+1)*sizeof(OBJECT));
334     else
335     rptr = r_set;
336     if(!rptr)
337     goto memerr;
338     setintersect(rptr,del_set,optr);
339    
340     if(QT_SET_CNT(rptr) > 0)
341     {
342     /* First need to check if set size larger than QT_MAXSET: if so
343     need to allocate larger array
344     */
345     sptr = QT_SET_PTR(rptr);
346     for(i = QT_SET_CNT(rptr); i > 0; i--)
347     {
348     t_id = QT_SET_NEXT_ELEM(sptr);
349     qt = qtdelelem(qt,t_id);
350     }
351     }
352     /* If we allocated memory: free it */
353     if(rptr != r_set)
354     free(rptr);
355    
356     return(qt);
357     memerr:
358     error(SYSTEM,"delete_tri_set():Unable to allocate memory");
359     }
360    
361     QUADTREE
362     expand_node(qt,q0,q1,q2,optr,n)
363     QUADTREE qt;
364     FVECT q0,q1,q2;
365     OBJECT *optr;
366     int n;
367     {
368     OBJECT *tptr,t_set[QT_MAXSET+1];
369     int i,t_id,found;
370     TRI *t;
371     FVECT v0,v1,v2;
372    
373     if(QT_SET_CNT(optr) > QT_MAXSET)
374     tptr = (OBJECT *)malloc((QT_SET_CNT(optr)+1)*sizeof(OBJECT));
375     else
376     tptr = t_set;
377     if(!tptr)
378     goto memerr;
379    
380     qtgetset(tptr,qt);
381     /* If set size exceeds threshold: subdivide cell and reinsert tris*/
382     qtfreeleaf(qt);
383     qtSubdivide(qt);
384    
385     for(optr = QT_SET_PTR(tptr),i=QT_SET_CNT(tptr); i > 0; i--)
386     {
387     t_id = QT_SET_NEXT_ELEM(optr);
388     t = SM_NTH_TRI(smMesh,t_id);
389     if(!T_IS_VALID(t))
390     continue;
391     VSUB(v0,SM_T_NTH_WV(smMesh,t,0),SM_VIEW_CENTER(smMesh));
392     VSUB(v1,SM_T_NTH_WV(smMesh,t,1),SM_VIEW_CENTER(smMesh));
393     VSUB(v2,SM_T_NTH_WV(smMesh,t,2),SM_VIEW_CENTER(smMesh));
394     qt = qtAdd_tri(qt,q0,q1,q2,v0,v1,v2,t_id,n);
395     }
396     /* If we allocated memory: free it */
397     if( tptr != t_set)
398     free(tptr);
399    
400     return(qt);
401     memerr:
402     error(SYSTEM,"expand_node():Unable to allocate memory");
403     }
404    
405     add_tri_expand(qtptr,f,argptr,q0,q1,q2,t0,t1,t2,n)
406 gwlarson 3.5 QUADTREE *qtptr;
407 gwlarson 3.8 int *f;
408     ADD_ARGS *argptr;
409 gwlarson 3.5 FVECT q0,q1,q2;
410     FVECT t0,t1,t2;
411     int n;
412     {
413 gwlarson 3.8 int t_id;
414     OBJECT *optr,*del_set;
415 gwlarson 3.7
416 gwlarson 3.8 t_id = argptr->t_id;
417    
418     if(QT_IS_EMPTY(*qtptr))
419     {
420     *qtptr = qtaddelem(*qtptr,t_id);
421     return;
422     }
423     if(!QT_LEAF_IS_FLAG(*qtptr))
424     {
425 gwlarson 3.5 optr = qtqueryset(*qtptr);
426    
427 gwlarson 3.8 if(del_set=argptr->del_set)
428     *qtptr = delete_tri_set(*qtptr,optr,del_set);
429     *qtptr = qtaddelem(*qtptr,t_id);
430     }
431     if (n >= QT_MAX_LEVELS)
432     return;
433     optr = qtqueryset(*qtptr);
434     if(QT_SET_CNT(optr) < QT_SET_THRESHOLD)
435     return;
436     *qtptr = expand_node(*qtptr,q0,q1,q2,optr,n);
437 gwlarson 3.5 }
438    
439    
440 gwlarson 3.8
441     add_tri(qtptr,fptr,argptr)
442 gwlarson 3.5 QUADTREE *qtptr;
443     int *fptr;
444 gwlarson 3.8 ADD_ARGS *argptr;
445 gwlarson 3.5 {
446    
447 gwlarson 3.8 OBJECT *optr,*del_set;
448     int t_id;
449    
450     t_id = argptr->t_id;
451    
452    
453     if(QT_IS_EMPTY(*qtptr))
454     {
455     *qtptr = qtaddelem(*qtptr,t_id);
456     if(!QT_FLAG_FILL_TRI(*fptr))
457     (*fptr)++;
458     return;
459     }
460     if(QT_LEAF_IS_FLAG(*qtptr))
461     return;
462    
463     optr = qtqueryset(*qtptr);
464    
465     if(del_set = argptr->del_set)
466     *qtptr = delete_tri_set(*qtptr,optr,del_set);
467    
468     if(!QT_IS_EMPTY(*qtptr))
469 gwlarson 3.5 {
470 gwlarson 3.8 optr = qtqueryset(*qtptr);
471     if(QT_SET_CNT(optr) >= QT_SET_THRESHOLD)
472     (*fptr) |= QT_EXPAND;
473 gwlarson 3.5 }
474 gwlarson 3.8 if(!QT_FLAG_FILL_TRI(*fptr))
475     (*fptr)++;
476     *qtptr = qtaddelem(*qtptr,t_id);
477    
478 gwlarson 3.5 }
479    
480    
481 gwlarson 3.7 smLocator_add_tri(sm,t_id,v0_id,v1_id,v2_id,del_set)
482 gwlarson 3.1 SM *sm;
483     int t_id;
484 gwlarson 3.7 int v0_id,v1_id,v2_id;
485     OBJECT *del_set;
486 gwlarson 3.1 {
487     STREE *st;
488 gwlarson 3.5 FVECT v0,v1,v2;
489 gwlarson 3.8 ADD_ARGS args;
490    
491 gwlarson 3.1 st = SM_LOCATOR(sm);
492    
493 gwlarson 3.8
494 gwlarson 3.5 VSUB(v0,SM_NTH_WV(sm,v0_id),SM_VIEW_CENTER(sm));
495     VSUB(v1,SM_NTH_WV(sm,v1_id),SM_VIEW_CENTER(sm));
496     VSUB(v2,SM_NTH_WV(sm,v2_id),SM_VIEW_CENTER(sm));
497 gwlarson 3.1
498 gwlarson 3.7 qtClearAllFlags();
499 gwlarson 3.8 args.t_id = t_id;
500     args.del_set = del_set;
501 gwlarson 3.1
502 gwlarson 3.8 stApply_to_tri(st,v0,v1,v2,add_tri,add_tri_expand,&args);
503    
504 gwlarson 3.1 }
505    
506     /* Add a triangle to the base array with vertices v1-v2-v3 */
507     int
508 gwlarson 3.8 smAdd_tri(sm, v0_id,v1_id,v2_id)
509 gwlarson 3.1 SM *sm;
510     int v0_id,v1_id,v2_id;
511     {
512     int t_id;
513     TRI *t;
514    
515 gwlarson 3.8 t_id = smAlloc_tri(sm);
516 gwlarson 3.1
517     if(t_id == -1)
518     return(t_id);
519    
520     t = SM_NTH_TRI(sm,t_id);
521 gwlarson 3.8
522 gwlarson 3.1 T_CLEAR_NBRS(t);
523 gwlarson 3.8 /* set the triangle vertex ids */
524     T_NTH_V(t,0) = v0_id;
525     T_NTH_V(t,1) = v1_id;
526     T_NTH_V(t,2) = v2_id;
527 gwlarson 3.1
528 gwlarson 3.8 SM_NTH_VERT(sm,v0_id) = t_id;
529     SM_NTH_VERT(sm,v1_id) = t_id;
530     SM_NTH_VERT(sm,v2_id) = t_id;
531    
532 gwlarson 3.1 if(SM_BASE_ID(sm,v0_id) || SM_BASE_ID(sm,v1_id) || SM_BASE_ID(sm,v2_id))
533     {
534     smClear_tri_flags(sm,t_id);
535     SM_SET_NTH_T_BASE(sm,t_id);
536     }
537     else
538     {
539 gwlarson 3.8 SM_CLR_NTH_T_BASE(sm,t_id);
540 gwlarson 3.1 SM_SET_NTH_T_ACTIVE(sm,t_id);
541     SM_SET_NTH_T_NEW(sm,t_id);
542 gwlarson 3.8 S_SET_FLAG(T_NTH_V(t,0));
543     S_SET_FLAG(T_NTH_V(t,1));
544     S_SET_FLAG(T_NTH_V(t,2));
545     SM_SAMPLE_TRIS(sm)++;
546 gwlarson 3.1 smNew_tri_cnt++;
547     }
548    
549     /* return initialized triangle */
550     return(t_id);
551     }
552    
553    
554     void
555 gwlarson 3.8 smTris_swap_edge(sm,t_id,t1_id,e,e1,tn_id,tn1_id,add_ptr,delptr)
556 gwlarson 3.1 SM *sm;
557     int t_id,t1_id;
558     int e,e1;
559 gwlarson 3.5 int *tn_id,*tn1_id;
560 gwlarson 3.7 LIST **add_ptr;
561 gwlarson 3.8 QUADTREE *delptr;
562 gwlarson 3.1
563     {
564     int verts[3],enext,eprev,e1next,e1prev;
565     TRI *n;
566     FVECT p1,p2,p3;
567     int ta_id,tb_id;
568     /* swap diagonal (e relative to t, and e1 relative to t1)
569     defined by quadrilateral
570     formed by t,t1- swap for the opposite diagonal
571     */
572     enext = (e+1)%3;
573     eprev = (e+2)%3;
574     e1next = (e1+1)%3;
575     e1prev = (e1+2)%3;
576 gwlarson 3.8 verts[e] = T_NTH_V(SM_NTH_TRI(sm,t_id),e);
577     verts[enext] = T_NTH_V(SM_NTH_TRI(sm,t1_id),e1prev);
578     verts[eprev] = T_NTH_V(SM_NTH_TRI(sm,t_id),eprev);
579     ta_id = smAdd_tri(sm,verts[0],verts[1],verts[2]);
580 gwlarson 3.7 *add_ptr = push_data(*add_ptr,ta_id);
581 gwlarson 3.8 verts[e1] = T_NTH_V(SM_NTH_TRI(sm,t1_id),e1);
582     verts[e1next] = T_NTH_V(SM_NTH_TRI(sm,t_id),eprev);
583     verts[e1prev] = T_NTH_V(SM_NTH_TRI(sm,t1_id),e1prev);
584     tb_id = smAdd_tri(sm,verts[0],verts[1],verts[2]);
585 gwlarson 3.7 *add_ptr = push_data(*add_ptr,tb_id);
586 gwlarson 3.3
587 gwlarson 3.1 /* set the neighbors */
588 gwlarson 3.8 T_NTH_NBR(SM_NTH_TRI(sm,ta_id),e) = T_NTH_NBR(SM_NTH_TRI(sm,t1_id),e1next);
589     T_NTH_NBR(SM_NTH_TRI(sm,tb_id),e1) = T_NTH_NBR(SM_NTH_TRI(sm,t_id),enext);
590     T_NTH_NBR(SM_NTH_TRI(sm,ta_id),enext) =tb_id;
591     T_NTH_NBR(SM_NTH_TRI(sm,tb_id),e1next) = ta_id;
592     T_NTH_NBR(SM_NTH_TRI(sm,ta_id),eprev)=T_NTH_NBR(SM_NTH_TRI(sm,t_id),eprev);
593     T_NTH_NBR(SM_NTH_TRI(sm,tb_id),e1prev)=
594     T_NTH_NBR(SM_NTH_TRI(sm,t1_id),e1prev);
595 gwlarson 3.1
596     /* Reset neighbor pointers of original neighbors */
597 gwlarson 3.8 n = SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,t_id),enext));
598 gwlarson 3.1 T_NTH_NBR(n,T_NTH_NBR_PTR(t_id,n)) = tb_id;
599 gwlarson 3.8 n = SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,t_id),eprev));
600 gwlarson 3.1 T_NTH_NBR(n,T_NTH_NBR_PTR(t_id,n)) = ta_id;
601    
602 gwlarson 3.8 n = SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,t1_id),e1next));
603 gwlarson 3.1 T_NTH_NBR(n,T_NTH_NBR_PTR(t1_id,n)) = ta_id;
604 gwlarson 3.8 n = SM_NTH_TRI(sm,T_NTH_NBR(SM_NTH_TRI(sm,t1_id),e1prev));
605 gwlarson 3.1 T_NTH_NBR(n,T_NTH_NBR_PTR(t1_id,n)) = tb_id;
606    
607     /* Delete two parent triangles */
608 gwlarson 3.7 if(remove_from_list(t_id,add_ptr))
609     smDelete_tri(sm,t_id);
610     else
611 gwlarson 3.8 *delptr = qtaddelem(*delptr,t_id);
612 gwlarson 3.5
613 gwlarson 3.7 if(remove_from_list(t1_id,add_ptr))
614     smDelete_tri(sm,t1_id);
615     else
616 gwlarson 3.8 *delptr = qtaddelem(*delptr,t1_id);
617 gwlarson 3.7
618 gwlarson 3.1 *tn_id = ta_id;
619     *tn1_id = tb_id;
620     }
621    
622 gwlarson 3.7 smUpdate_locator(sm,add_list,del_set)
623 gwlarson 3.3 SM *sm;
624 gwlarson 3.7 LIST *add_list;
625     OBJECT *del_set;
626 gwlarson 3.3 {
627 gwlarson 3.7 int t_id,i;
628 gwlarson 3.3 TRI *t;
629 gwlarson 3.7 OBJECT *optr;
630    
631 gwlarson 3.3 while(add_list)
632     {
633     t_id = pop_list(&add_list);
634 gwlarson 3.6 t = SM_NTH_TRI(sm,t_id);
635 gwlarson 3.7 smLocator_add_tri(sm,t_id,T_NTH_V(t,0),T_NTH_V(t,1),T_NTH_V(t,2),del_set);
636 gwlarson 3.3 }
637 gwlarson 3.7
638     optr = QT_SET_PTR(del_set);
639     for(i = QT_SET_CNT(del_set); i > 0; i--)
640 gwlarson 3.3 {
641 gwlarson 3.7 t_id = QT_SET_NEXT_ELEM(optr);
642 gwlarson 3.8 #if 0
643     t = SM_NTH_TRI(sm,t_id);
644     smLocator_remove_tri(sm,t_id,T_NTH_V(t,0),T_NTH_V(t,1),T_NTH_V(t,2));
645     #endif
646 gwlarson 3.7 smDelete_tri(sm,t_id);
647 gwlarson 3.3 }
648     }
649 gwlarson 3.1 /* MUST add check for constrained edges */
650 gwlarson 3.3 int
651 gwlarson 3.8 smFix_tris(sm,id,tlist,add_list,delptr)
652 gwlarson 3.1 SM *sm;
653     int id;
654 gwlarson 3.8 LIST *tlist,*add_list;
655     QUADTREE *delptr;
656 gwlarson 3.1 {
657     TRI *t,*t_opp;
658 gwlarson 3.3 FVECT p,p1,p2,p3;
659     int e,e1,swapped = 0;
660 gwlarson 3.1 int t_id,t_opp_id;
661 gwlarson 3.3
662     VSUB(p,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
663 gwlarson 3.1 while(tlist)
664     {
665 gwlarson 3.3 t_id = pop_list(&tlist);
666 gwlarson 3.1 t = SM_NTH_TRI(sm,t_id);
667     e = (T_WHICH_V(t,id)+1)%3;
668     t_opp_id = T_NTH_NBR(t,e);
669     t_opp = SM_NTH_TRI(sm,t_opp_id);
670 gwlarson 3.8
671 gwlarson 3.7 smDir_in_cone(sm,p1,T_NTH_V(t_opp,0));
672     smDir_in_cone(sm,p2,T_NTH_V(t_opp,1));
673     smDir_in_cone(sm,p3,T_NTH_V(t_opp,2));
674 gwlarson 3.1 if(point_in_cone(p,p1,p2,p3))
675     {
676     swapped = 1;
677     e1 = T_NTH_NBR_PTR(t_id,t_opp);
678     /* check list for t_opp and Remove if there */
679     remove_from_list(t_opp_id,&tlist);
680 gwlarson 3.3 smTris_swap_edge(sm,t_id,t_opp_id,e,e1,&t_id,&t_opp_id,
681 gwlarson 3.8 &add_list,delptr);
682 gwlarson 3.1 tlist = push_data(tlist,t_id);
683     tlist = push_data(tlist,t_opp_id);
684     }
685     }
686 gwlarson 3.8 smUpdate_locator(sm,add_list,qtqueryset(*delptr));
687 gwlarson 3.1 return(swapped);
688     }
689    
690     /* Give the vertex "id" and a triangle "t" that it belongs to- return the
691     next nbr in a counter clockwise order about vertex "id"
692     */
693     int
694     smTri_next_ccw_nbr(sm,t,id)
695     SM *sm;
696     TRI *t;
697     int id;
698     {
699     int t_id;
700 gwlarson 3.8 int nbr_id;
701 gwlarson 3.1
702     /* Want the edge for which "id" is the destination */
703     t_id = (T_WHICH_V(t,id)+ 2)% 3;
704 gwlarson 3.8 nbr_id = T_NTH_NBR(t,t_id);
705     return(nbr_id);
706 gwlarson 3.1 }
707    
708     smClear_tri_flags(sm,id)
709     SM *sm;
710     int id;
711     {
712     int i;
713    
714     for(i=0; i < T_FLAGS; i++)
715 gwlarson 3.8 SM_CLR_NTH_T_FLAG(sm,id,i);
716 gwlarson 3.1
717     }
718    
719     /* Locate the point-id in the point location structure: */
720     int
721 gwlarson 3.8 smInsert_samp(sm,s_id,tri_id)
722 gwlarson 3.1 SM *sm;
723     int s_id,tri_id;
724     {
725 gwlarson 3.8 int v0_id,v1_id,v2_id;
726     int t0_id,t1_id,t2_id,replaced;
727 gwlarson 3.7 LIST *tlist,*add_list;
728 gwlarson 3.8 OBJECT del_set[2];
729     QUADTREE delnode;
730 gwlarson 3.1 FVECT npt;
731 gwlarson 3.8 TRI *tri,*nbr;
732 gwlarson 3.1
733 gwlarson 3.7 add_list = NULL;
734 gwlarson 3.1
735 gwlarson 3.8 v0_id = T_NTH_V(SM_NTH_TRI(sm,tri_id),0);
736     v1_id = T_NTH_V(SM_NTH_TRI(sm,tri_id),1);
737     v2_id = T_NTH_V(SM_NTH_TRI(sm,tri_id),2);
738    
739     t0_id = smAdd_tri(sm,s_id,v0_id,v1_id);
740 gwlarson 3.3 /* Add triangle to the locator */
741 gwlarson 3.6
742     add_list = push_data(add_list,t0_id);
743 gwlarson 3.3
744 gwlarson 3.8 t1_id = smAdd_tri(sm,s_id,v1_id,v2_id);
745 gwlarson 3.6 add_list = push_data(add_list,t1_id);
746    
747 gwlarson 3.8 t2_id = smAdd_tri(sm,s_id,v2_id,v0_id);
748 gwlarson 3.6 add_list = push_data(add_list,t2_id);
749 gwlarson 3.3
750 gwlarson 3.8 /* CAUTION: tri must be assigned after Add_tri: pointers may change*/
751     tri = SM_NTH_TRI(sm,tri_id);
752    
753 gwlarson 3.1 /* Set the neighbor pointers for the new tris */
754 gwlarson 3.8 T_NTH_NBR(SM_NTH_TRI(sm,t0_id),0) = t2_id;
755     T_NTH_NBR(SM_NTH_TRI(sm,t0_id),1) = T_NTH_NBR(tri,0);
756     T_NTH_NBR(SM_NTH_TRI(sm,t0_id),2) = t1_id;
757     T_NTH_NBR(SM_NTH_TRI(sm,t1_id),0) = t0_id;
758     T_NTH_NBR(SM_NTH_TRI(sm,t1_id),1) = T_NTH_NBR(tri,1);
759     T_NTH_NBR(SM_NTH_TRI(sm,t1_id),2) = t2_id;
760     T_NTH_NBR(SM_NTH_TRI(sm,t2_id),0) = t1_id;
761     T_NTH_NBR(SM_NTH_TRI(sm,t2_id),1) = T_NTH_NBR(tri,2);
762     T_NTH_NBR(SM_NTH_TRI(sm,t2_id),2) = t0_id;
763 gwlarson 3.1
764     /* Reset the neigbor pointers for the neighbors of the original */
765     nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,0));
766     T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t0_id;
767     nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,1));
768     T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t1_id;
769     nbr = SM_NTH_TRI(sm,T_NTH_NBR(tri,2));
770     T_NTH_NBR(nbr,T_NTH_NBR_PTR(tri_id,nbr)) = t2_id;
771    
772 gwlarson 3.8 del_set[0] = 1; del_set[1] = tri_id;
773     delnode = qtnewleaf(del_set);
774 gwlarson 3.6
775 gwlarson 3.1 /* Fix up the new triangles*/
776     tlist = push_data(NULL,t0_id);
777     tlist = push_data(tlist,t1_id);
778     tlist = push_data(tlist,t2_id);
779    
780 gwlarson 3.8 smFix_tris(sm,s_id,tlist,add_list,&delnode);
781 gwlarson 3.1
782 gwlarson 3.8 qtfreeleaf(delnode);
783 gwlarson 3.1
784 gwlarson 3.8 return(TRUE);
785 gwlarson 3.1 }
786    
787    
788     int
789 gwlarson 3.8 smTri_in_set(sm,p,optr)
790 gwlarson 3.1 SM *sm;
791 gwlarson 3.8 FVECT p;
792     OBJECT *optr;
793 gwlarson 3.1 {
794 gwlarson 3.8 int i,t_id;
795     FVECT n,v0,v1,v2;
796     TRI *t;
797 gwlarson 3.1
798 gwlarson 3.8 for (i = QT_SET_CNT(optr),optr = QT_SET_PTR(optr);i > 0; i--)
799 gwlarson 3.1 {
800 gwlarson 3.8 /* Find the first triangle that pt falls */
801     t_id = QT_SET_NEXT_ELEM(optr);
802     t = SM_NTH_TRI(sm,t_id);
803     if(!T_IS_VALID(t))
804     continue;
805     VSUB(v0,SM_T_NTH_WV(sm,t,0),SM_VIEW_CENTER(sm));
806     VSUB(v1,SM_T_NTH_WV(sm,t,1),SM_VIEW_CENTER(sm));
807     VSUB(v2,SM_T_NTH_WV(sm,t,2),SM_VIEW_CENTER(sm));
808    
809     if(EQUAL_VEC3(v0,p) || EQUAL_VEC3(v1,p) || EQUAL_VEC3(v2,p))
810     return(t_id);
811    
812     VCROSS(n,v1,v0);
813     if(DOT(n,p) >0.0)
814     continue;
815     VCROSS(n,v2,v1);
816     if(DOT(n,p) > 0.0)
817     continue;
818    
819     VCROSS(n,v0,v2);
820     if(DOT(n,p) > 0.0)
821     continue;
822    
823     return(t_id);
824 gwlarson 3.1 }
825 gwlarson 3.8 return(INVALID);
826 gwlarson 3.1 }
827    
828 gwlarson 3.8 int
829     smPointLocateTri(sm,id)
830 gwlarson 3.1 SM *sm;
831 gwlarson 3.8 int id;
832 gwlarson 3.1 {
833 gwlarson 3.8 FVECT tpt;
834     QUADTREE qt,*optr;
835    
836     VSUB(tpt,SM_NTH_WV(sm,id),SM_VIEW_CENTER(sm));
837     qt = smPointLocateCell(sm,tpt);
838     if(QT_IS_EMPTY(qt))
839     return(INVALID);
840    
841     optr = qtqueryset(qt);
842     return(smTri_in_set(sm,tpt,optr));
843     }
844 gwlarson 3.1
845 gwlarson 3.8
846     /*
847     Determine whether this sample should:
848     a) be added to the mesh by subdividing the triangle
849     b) ignored
850     c) replace values of an existing mesh vertex
851    
852     In case a, the routine will return TRUE, for b,c will return FALSE
853     In case a, also determine if this sample should cause the deletion of
854     an existing vertex. If so *rptr will contain the id of the sample to delete
855     after the new sample has been added.
856    
857     ASSUMPTION: this will not be called with a new sample that is
858     a BASE point.
859    
860     The following tests are performed (in order) to determine the fate
861     of the sample:
862    
863     1) If the world space point of the sample coincides with one of the
864     triangle vertex samples: The values of the triangle sample are
865     replaced with the new values and FALSE is returned
866     2) If the new sample is close in ws, and close in the spherical projection
867     to one of the triangle vertex samples:
868     pick the point with dir closest to that of the canonical view
869     If it is the new point: mark existing point for deletion,and return
870     TRUE,else return FALSE
871     3) If the spherical projection of new is < REPLACE_EPS from a base point:
872     add new and mark the base for deletion: return TRUE
873     4) If the addition of the new sample point would introduce a "puncture"
874     or cause new triangles with large depth differences:return FALSE
875     This attempts to throw out points that should be occluded
876     */
877     int
878     smTest_sample(sm,tri_id,s_id,dir,rptr)
879     SM *sm;
880     int tri_id,s_id;
881     FVECT dir;
882     int *rptr;
883     {
884     TRI *tri;
885     double d,d2,dnear,dspt,d01,d12,d20,s0,s1,s2,ds,dv;
886     int vid[3],i,nearid,norm,dcnt,bcnt;
887     FVECT diff[3],spt,npt;
888     FVECT p;
889    
890     VCOPY(p,SM_NTH_WV(sm,s_id));
891     *rptr = INVALID;
892     bcnt = dcnt = 0;
893     tri = SM_NTH_TRI(sm,tri_id);
894     vid[0] = T_NTH_V(tri,0);
895     vid[1] = T_NTH_V(tri,1);
896     vid[2] = T_NTH_V(tri,2);
897    
898     /* TEST 1: Test if the new point has the same world space point as
899     one of the existing triangle vertices
900     */
901     for(i=0; i<3; i++)
902     {
903     if(SM_BASE_ID(sm,vid[i]))
904     {
905     bcnt++;
906     continue;
907     }
908     if(SM_DIR_ID(sm,vid[i]))
909     dcnt++;
910     VSUB(diff[i],SM_NTH_WV(sm,vid[i]),p);
911     /* If same world point: replace */
912     if(ZERO_VEC3(diff[i]))
913     {
914     sReset_samp(SM_SAMP(sm),vid[i],s_id);
915     SM_TONE_MAP(sm) = 0;
916     return(FALSE);
917     }
918     }
919     if(SM_DIR_ID(sm,s_id))
920     return(TRUE);
921     /* TEST 2: If the new sample is close in ws, and close in the spherical
922     projection to one of the triangle vertex samples
923     */
924     norm = FALSE;
925     if(bcnt + dcnt != 3)
926     {
927     VSUB(spt,p,SM_VIEW_CENTER(sm));
928     ds = DOT(spt,spt);
929     dnear = FHUGE;
930     for(i=0; i<3; i++)
931     {
932     if(SM_BASE_ID(sm,vid[i]) || SM_DIR_ID(sm,vid[i]))
933     continue;
934     d = DOT(diff[i],diff[i])/ds;
935     if(d < dnear)
936     {
937     dnear = d;
938     nearid=vid[i];
939     }
940     }
941    
942     if(dnear <= smMinSampDiff*smMinSampDiff)
943     {
944     /* Pick the point with dir closest to that of the canonical view
945     if it is the new sample: mark existing point for deletion
946     */
947     normalize(spt);
948     norm = TRUE;
949     VSUB(npt,SM_NTH_WV(sm,nearid),SM_VIEW_CENTER(sm));
950     normalize(npt);
951     d = fdir2diff(SM_NTH_W_DIR(sm,nearid), npt);
952     if(dir)
953     d2 = 2. - 2.*DOT(dir,spt);
954     else
955     d2 = fdir2diff(SM_NTH_W_DIR(sm,s_id), spt);
956     /* The existing sample is a better sample:punt */
957     if(d2 > d)
958     return(FALSE);
959     else
960     {
961     /* The new sample is better: mark the existing one
962     for deletion after the new one is added*/
963     *rptr = nearid;
964     return(TRUE);
965     }
966     }
967     }
968     /* TEST 3: If the spherical projection of new is < S_REPLACE_EPS
969     from a base point
970     */
971     if(bcnt)
972     {
973     dnear = FHUGE;
974     if(bcnt + dcnt ==3)
975     VSUB(spt,p,SM_VIEW_CENTER(sm));
976     if(!norm)
977     normalize(spt);
978    
979     for(i=0; i<3; i++)
980     {
981     if(!SM_BASE_ID(sm,vid[i]))
982     continue;
983     VSUB(npt,SM_NTH_WV(sm,vid[i]),SM_VIEW_CENTER(sm));
984     d = DIST_SQ(npt,spt);
985     if(d < S_REPLACE_EPS && d < dnear)
986     {
987     dnear = d;
988     nearid = vid[i];
989     }
990     }
991     if(dnear != FHUGE)
992     {
993     /* add new and mark the base for deletion */
994     *rptr = nearid;
995     return(TRUE);
996     }
997     }
998    
999     /* TEST 4:
1000     If the addition of the new sample point would introduce a "puncture"
1001     or cause new triangles with large depth differences:dont add:
1002     */
1003     if(bcnt || dcnt)
1004     return(TRUE);
1005     /* If the new point is in front of the existing points- add */
1006     dv = DIST_SQ(SM_NTH_WV(sm,vid[0]),SM_VIEW_CENTER(sm));
1007     if(ds < dv)
1008     return(TRUE);
1009    
1010     d01 = DIST_SQ(SM_NTH_WV(sm,vid[1]),SM_NTH_WV(sm,vid[0]));
1011     s0 = DOT(diff[0],diff[0]);
1012     if(s0 < S_REPLACE_SCALE*d01)
1013     return(TRUE);
1014     d12 = DIST_SQ(SM_NTH_WV(sm,vid[2]),SM_NTH_WV(sm,vid[1]));
1015     if(s0 < S_REPLACE_SCALE*d12)
1016     return(TRUE);
1017     d20 = DIST_SQ(SM_NTH_WV(sm,vid[0]),SM_NTH_WV(sm,vid[2]));
1018     if(s0 < S_REPLACE_SCALE*d20)
1019     return(TRUE);
1020     d = MIN3(d01,d12,d20);
1021     s1 = DOT(diff[1],diff[1]);
1022     if(s1 < S_REPLACE_SCALE*d)
1023     return(TRUE);
1024     s2 = DOT(diff[2],diff[2]);
1025     if(s2 < S_REPLACE_SCALE*d)
1026     return(TRUE);
1027    
1028    
1029     return(FALSE);
1030     }
1031    
1032    
1033     int
1034     smAlloc_samp(sm,c,dir,pt)
1035     SM *sm;
1036     COLR c;
1037     FVECT dir,pt;
1038     {
1039     int s_id,replaced,cnt;
1040     SAMP *s;
1041     FVECT p;
1042    
1043     s = SM_SAMP(sm);
1044     s_id = sAlloc_samp(s,&replaced);
1045    
1046     cnt=0;
1047     while(replaced)
1048 gwlarson 3.1 {
1049 gwlarson 3.8 if(smMesh_remove_vertex(sm,s_id))
1050     break;
1051     s_id = sAlloc_samp(s,&replaced);
1052     cnt++;
1053     if(cnt > S_MAX_SAMP(s))
1054     error(CONSISTENCY,"smAlloc_samp():unable to find free samp\n");
1055     }
1056 gwlarson 3.1
1057 gwlarson 3.8 if(pt)
1058     sInit_samp(s,s_id,c,dir,pt);
1059 gwlarson 3.1 else
1060 gwlarson 3.8 {
1061     VADD(p,dir,SM_VIEW_CENTER(sm));
1062     sInit_samp(s,s_id,c,NULL,p);
1063     }
1064     return(s_id);
1065 gwlarson 3.1 }
1066    
1067     int
1068 gwlarson 3.8 smAdd_samp(sm,s_id,p,dir)
1069 gwlarson 3.1 SM *sm;
1070     int s_id;
1071 gwlarson 3.8 FVECT p,dir;
1072 gwlarson 3.1 {
1073 gwlarson 3.8 int t_id,r_id,test;
1074 gwlarson 3.1 double d;
1075 gwlarson 3.8
1076     r_id = INVALID;
1077     t_id = smPointLocateTri(sm,s_id);
1078     if(t_id == INVALID)
1079 gwlarson 3.1 {
1080     #ifdef DEBUG
1081 gwlarson 3.8 eputs("smAdd_samp(): tri not found \n");
1082 gwlarson 3.1 #endif
1083 gwlarson 3.8 return(FALSE);
1084 gwlarson 3.1 }
1085 gwlarson 3.8 if(!SM_BASE_ID(sm,s_id))
1086 gwlarson 3.1 {
1087 gwlarson 3.8 if(!smTest_sample(sm,t_id,s_id,dir,&r_id))
1088     return(FALSE);
1089     /* If not a sky point, add distance from the viewcenter to average*/
1090     if( !SM_DIR_ID(sm,s_id))
1091     {
1092     d = DIST(p,SM_VIEW_CENTER(smMesh));
1093     smDist_sum += 1.0/d;
1094     }
1095 gwlarson 3.1 }
1096 gwlarson 3.8 test = smInsert_samp(smMesh,s_id,t_id);
1097 gwlarson 3.5
1098 gwlarson 3.8 if(test && r_id != INVALID)
1099     smMesh_remove_vertex(sm,r_id);
1100    
1101     return(test);
1102 gwlarson 3.1 }
1103    
1104     /*
1105     * int
1106     * smNewSamp(c, dir, p) : register new sample point and return index
1107     * COLR c; : pixel color (RGBE)
1108     * FVECT dir; : ray direction vector
1109     * FVECT p; : world intersection point
1110     *
1111     * Add new sample point to data structures, removing old values as necessary.
1112     * New sample representation will be output in next call to smUpdate().
1113     * If the point is a sky point: then v=NULL
1114     */
1115     int
1116     smNewSamp(c,dir,p)
1117     COLR c;
1118     FVECT dir;
1119     FVECT p;
1120    
1121     {
1122     int s_id;
1123 gwlarson 3.8 int debug=0;
1124     static FILE *fp;
1125     static int cnt=0,n=3010;
1126 gwlarson 3.1 /* First check if this the first sample: if so initialize mesh */
1127     if(SM_NUM_SAMP(smMesh) == 0)
1128 gwlarson 3.8 {
1129     smInit_sm(smMesh,odev.v.vp);
1130     #if 0
1131     fp = fopen("Debug_data.view","w");
1132     fprintf(fp,"%d %d %f %f %f ",1280,1024,odev.v.vp[0],odev.v.vp[1],
1133     odev.v.vp[2]);
1134     fprintf(fp,"%f %f %f ",odev.v.vdir[0],odev.v.vdir[1],
1135     odev.v.vdir[2]);
1136     fprintf(fp,"%f %f %f ",odev.v.vup[0],odev.v.vup[1],odev.v.vup[2]);
1137     fprintf(fp,"%f %f ",odev.v.horiz,odev.v.vert);
1138     fclose(fp);
1139     fp = fopen("Debug_data","w");
1140     #endif
1141     }
1142     #if 0
1143     fprintf(fp,"%f %f %f %f %f %f ",p[0],p[1],p[2],(float)c[0]/255.0,(float)c[1]/255.0,
1144     (float)c[2]/255.0);
1145     #endif
1146    
1147     /* Allocate space for a sample */
1148     s_id = smAlloc_samp(smMesh,c,dir,p);
1149     #if 0
1150     if(cnt==n)
1151     return(-1);
1152     cnt++;
1153     #endif
1154     /* Add the sample to the mesh */
1155     if(!smAdd_samp(smMesh,s_id,p,dir))
1156     {
1157     /* If the sample space was not used: return */
1158     smUnalloc_samp(smMesh,s_id);
1159     s_id = INVALID;
1160     }
1161 gwlarson 3.1 return(s_id);
1162    
1163     }
1164     int
1165 gwlarson 3.8 smAdd_base_vertex(sm,v)
1166 gwlarson 3.1 SM *sm;
1167 gwlarson 3.8 FVECT v;
1168 gwlarson 3.1 {
1169     int id;
1170    
1171     /* First add coordinate to the sample array */
1172 gwlarson 3.8 id = sAdd_base_point(SM_SAMP(sm),v);
1173     if(id == INVALID)
1174     return(INVALID);
1175 gwlarson 3.1 /* Initialize triangle pointer to -1 */
1176     smClear_vert(sm,id);
1177     return(id);
1178     }
1179    
1180    
1181    
1182     /* Initialize a the point location DAG based on a 6 triangle tesselation
1183     of the unit sphere centered on the view center. The DAG structure
1184     contains 6 roots: one for each initial base triangle
1185     */
1186    
1187     int
1188     smCreate_base_mesh(sm,type)
1189     SM *sm;
1190     int type;
1191     {
1192 gwlarson 3.8 int i,s_id,tri_id,nbr_id;
1193 gwlarson 3.1 int p[4],ids[4];
1194     int v0_id,v1_id,v2_id;
1195 gwlarson 3.8 FVECT d,pt,cntr,v0,v1,v2;
1196 gwlarson 3.1
1197     /* First insert the base vertices into the sample point array */
1198    
1199     for(i=0; i < 4; i++)
1200     {
1201 gwlarson 3.8 VCOPY(cntr,smDefault_base[i]);
1202 gwlarson 3.5 cntr[0] += .01;
1203     cntr[1] += .02;
1204     cntr[2] += .03;
1205 gwlarson 3.6 VADD(cntr,cntr,SM_VIEW_CENTER(sm));
1206 gwlarson 3.8 p[i] = smAdd_base_vertex(sm,cntr);
1207 gwlarson 3.1 }
1208     /* Create the base triangles */
1209     for(i=0;i < 4; i++)
1210     {
1211 gwlarson 3.8 v0_id = p[smTri_verts[i][0]];
1212     v1_id = p[smTri_verts[i][1]];
1213     v2_id = p[smTri_verts[i][2]];
1214     ids[i] = smAdd_tri(sm, v0_id,v1_id,v2_id);
1215 gwlarson 3.7 smLocator_add_tri(sm,ids[i],v0_id,v1_id,v2_id,NULL);
1216 gwlarson 3.1 }
1217 gwlarson 3.8
1218 gwlarson 3.1 /* Set neighbors */
1219    
1220 gwlarson 3.5 for(tri_id=0;tri_id < 4; tri_id++)
1221 gwlarson 3.8 for(nbr_id=0; nbr_id < 3; nbr_id++)
1222     T_NTH_NBR(SM_NTH_TRI(sm,ids[tri_id]),nbr_id) = smBase_nbrs[tri_id][nbr_id];
1223 gwlarson 3.1
1224 gwlarson 3.8
1225     /* Now add the centroids of the faces */
1226     for(tri_id=0;tri_id < 4; tri_id++)
1227     {
1228     VCOPY(v0,SM_T_NTH_WV(sm,SM_NTH_TRI(sm,ids[tri_id]),0));
1229     VCOPY(v1,SM_T_NTH_WV(sm,SM_NTH_TRI(sm,ids[tri_id]),1));
1230     VCOPY(v2,SM_T_NTH_WV(sm,SM_NTH_TRI(sm,ids[tri_id]),2));
1231     tri_centroid(v0,v1,v2,cntr);
1232     VSUB(cntr,cntr,SM_VIEW_CENTER(sm));
1233     normalize(cntr);
1234     VADD(cntr,cntr,SM_VIEW_CENTER(sm));
1235     s_id = smAdd_base_vertex(sm,cntr);
1236     smAdd_samp(sm,s_id,NULL,NULL);
1237     }
1238     return(1);
1239 gwlarson 3.1
1240     }
1241    
1242    
1243     int
1244     smNext_tri_flag_set(sm,i,which,b)
1245     SM *sm;
1246     int i,which;
1247 gwlarson 3.5 int b;
1248 gwlarson 3.1 {
1249    
1250 gwlarson 3.8 for(; i < SM_NUM_TRI(sm);i++)
1251 gwlarson 3.1 {
1252    
1253 gwlarson 3.8 if(!T_IS_VALID(SM_NTH_TRI(sm,i)))
1254 gwlarson 3.5 continue;
1255 gwlarson 3.8 if(!SM_IS_NTH_T_FLAG(sm,i,which))
1256     continue;
1257 gwlarson 3.1 if(!b)
1258     break;
1259     if((b==1) && !SM_BG_TRI(sm,i))
1260     break;
1261     if((b==2) && SM_BG_TRI(sm,i))
1262     break;
1263     }
1264    
1265     return(i);
1266     }
1267    
1268    
1269     int
1270     smNext_valid_tri(sm,i)
1271     SM *sm;
1272     int i;
1273     {
1274    
1275 gwlarson 3.8 while( i < SM_NUM_TRI(sm) && !T_IS_VALID(SM_NTH_TRI(sm,i)))
1276 gwlarson 3.1 i++;
1277    
1278     return(i);
1279     }
1280    
1281    
1282    
1283 gwlarson 3.8 qtTri_from_id(t_id,v0,v1,v2)
1284 gwlarson 3.1 int t_id;
1285     FVECT v0,v1,v2;
1286     {
1287     TRI *t;
1288    
1289     t = SM_NTH_TRI(smMesh,t_id);
1290 gwlarson 3.8 if(!T_IS_VALID(t))
1291     return(0);
1292     VSUB(v0,SM_T_NTH_WV(smMesh,t,0),SM_VIEW_CENTER(smMesh));
1293     VSUB(v1,SM_T_NTH_WV(smMesh,t,1),SM_VIEW_CENTER(smMesh));
1294     VSUB(v2,SM_T_NTH_WV(smMesh,t,2),SM_VIEW_CENTER(smMesh));
1295     return(1);
1296 gwlarson 3.1 }
1297    
1298    
1299 gwlarson 3.8 smRebuild_mesh(sm,v)
1300 gwlarson 3.1 SM *sm;
1301 gwlarson 3.8 VIEW *v;
1302 gwlarson 3.1 {
1303 gwlarson 3.8 int i,j,cnt;
1304     FVECT p,ov,dir;
1305     double d;
1306 gwlarson 3.1
1307 gwlarson 3.8 #ifdef DEBUG
1308     eputs("smRebuild_mesh(): rebuilding....");
1309     #endif
1310 gwlarson 3.1 /* Clear the mesh- and rebuild using the current sample array */
1311 gwlarson 3.8 /* Remember the current number of samples */
1312     cnt = SM_NUM_SAMP(sm);
1313     /* Calculate the difference between the current and new viewpoint*/
1314     /* Will need to subtract this off of sky points */
1315     VSUB(ov,v->vp,SM_VIEW_CENTER(sm));
1316     /* Initialize the mesh to 0 samples and the base triangles */
1317 gwlarson 3.1
1318 gwlarson 3.8 /* Go through all samples and add in if in the new view frustum, and
1319     the dir is <= 30 degrees from new view
1320     */
1321     j=0;
1322     for(i=0; i < cnt; i++)
1323 gwlarson 3.1 {
1324 gwlarson 3.8 /* First check if sample visible(conservative approx: if one of tris
1325     attached to sample is in frustum) */
1326     if(!S_IS_FLAG(i))
1327     continue;
1328    
1329     /* Next test if direction > 30 from current direction */
1330     if(SM_NTH_W_DIR(sm,i)!=-1)
1331     {
1332     VSUB(p,SM_NTH_WV(sm,i),v->vp);
1333     normalize(p);
1334     d = fdir2diff(SM_NTH_W_DIR(sm,i), p);
1335     if (d > MAXDIFF2)
1336     continue;
1337     }
1338     sReset_samp(SM_SAMP(sm),j,i);
1339     j++;
1340 gwlarson 3.1 }
1341 gwlarson 3.8 smInit_sm(sm,v->vp);
1342     for(i=0; i< j; i++)
1343     {
1344     S_SET_FLAG(i);
1345     VCOPY(p,SM_NTH_WV(sm,i));
1346     smAdd_samp(sm,i,p,NULL);
1347     }
1348     SM_NUM_SAMP(sm) = j;
1349     smNew_tri_cnt = SM_SAMPLE_TRIS(sm);
1350     #ifdef DEBUG
1351     eputs("smRebuild_mesh():done\n");
1352     #endif
1353 gwlarson 3.1 }
1354 gwlarson 3.5
1355     int
1356     intersect_tri_set(t_set,orig,dir,pt)
1357     OBJECT *t_set;
1358     FVECT orig,dir,pt;
1359     {
1360     OBJECT *optr;
1361     int i,t_id,id;
1362     int pid0,pid1,pid2;
1363     FVECT p0,p1,p2,p;
1364     TRI *t;
1365    
1366     optr = QT_SET_PTR(t_set);
1367     for(i = QT_SET_CNT(t_set); i > 0; i--)
1368     {
1369     t_id = QT_SET_NEXT_ELEM(optr);
1370    
1371     t = SM_NTH_TRI(smMesh,t_id);
1372 gwlarson 3.8 if(!T_IS_VALID(t))
1373     continue;
1374    
1375 gwlarson 3.5 pid0 = T_NTH_V(t,0);
1376     pid1 = T_NTH_V(t,1);
1377     pid2 = T_NTH_V(t,2);
1378     VCOPY(p0,SM_NTH_WV(smMesh,pid0));
1379     VCOPY(p1,SM_NTH_WV(smMesh,pid1));
1380     VCOPY(p2,SM_NTH_WV(smMesh,pid2));
1381     if(ray_intersect_tri(orig,dir,p0,p1,p2,p))
1382     {
1383     id = closest_point_in_tri(p0,p1,p2,p,pid0,pid1,pid2);
1384    
1385     if(pt)
1386     VCOPY(pt,p);
1387     #ifdef DEBUG_TEST_DRIVER
1388     Pick_tri = t_id;
1389     Pick_samp = id;
1390     VCOPY(Pick_point[0],p);
1391     #endif
1392     return(id);
1393     }
1394     }
1395     return(-1);
1396     }
1397    
1398 gwlarson 3.8 /* OS is constrained to be <= QT_MAXCSET : if the set exceeds this, the
1399     results of check_set are conservative
1400     */
1401    
1402     ray_trace_check_set(qtptr,fptr,argptr)
1403 gwlarson 3.5 QUADTREE *qtptr;
1404 gwlarson 3.8 int *fptr;
1405     RT_ARGS *argptr;
1406 gwlarson 3.5 {
1407 gwlarson 3.8 OBJECT tset[QT_MAXSET+1],*tptr;
1408 gwlarson 3.5 double dt,t;
1409     int found;
1410    
1411 gwlarson 3.8 if(QT_LEAF_IS_FLAG(*qtptr))
1412     {
1413     QT_FLAG_SET_DONE(*fptr);
1414     #if DEBUG
1415     eputs("ray_trace_check_set():Already visited this node:aborting\n");
1416     #endif
1417     return;
1418     }
1419 gwlarson 3.5 if(!QT_IS_EMPTY(*qtptr))
1420 gwlarson 3.8 {
1421     tptr = qtqueryset(*qtptr);
1422     if(QT_SET_CNT(tptr) > QT_MAXSET)
1423     tptr = (OBJECT *)malloc((QT_SET_CNT(tptr)+1)*sizeof(OBJECT));
1424     else
1425     tptr = tset;
1426     if(!tptr)
1427     goto memerr;
1428    
1429     qtgetset(tptr,*qtptr);
1430 gwlarson 3.5 /* Check triangles in set against those seen so far(os):only
1431     check new triangles for intersection (t_set')
1432     */
1433 gwlarson 3.8 check_set_large(tptr,argptr->os);
1434     if((found = intersect_tri_set(tptr,argptr->orig,argptr->dir,NULL))!= -1)
1435     {
1436     argptr->t_id = found;
1437     if(tptr != tset)
1438     free(tptr);
1439     QT_FLAG_SET_DONE(*fptr);
1440     return;
1441 gwlarson 3.5 }
1442 gwlarson 3.8 if(tptr != tset)
1443     free(tptr);
1444     }
1445     return;
1446     memerr:
1447     error(SYSTEM,"ray_trace_check_set():Unable to allocate memory");
1448 gwlarson 3.5 }
1449    
1450 gwlarson 3.8
1451     /*
1452     * int
1453     * smFindSamp(FVECT orig, FVECT dir)
1454     *
1455     * Find the closest sample to the given ray. Returns sample id, -1 on failure.
1456     * "dir" is assumed to be normalized
1457     */
1458    
1459 gwlarson 3.5 int
1460 gwlarson 3.6 smFindSamp(orig,dir)
1461 gwlarson 3.5 FVECT orig,dir;
1462     {
1463     FVECT b,p,o;
1464     OBJECT *ts;
1465     QUADTREE qt;
1466     int s_id;
1467     double d;
1468    
1469     /* r is the normalized vector from the view center to the current
1470     * ray point ( starting with "orig"). Find the cell that r falls in,
1471     * and test the ray against all triangles stored in the cell. If
1472     * the test fails, trace the projection of the ray across to the
1473     * next cell it intersects: iterate until either an intersection
1474     * is found, or the projection ray is // to the direction. The sample
1475     * corresponding to the triangle vertex closest to the intersection
1476     * point is returned.
1477     */
1478    
1479     /* First test if "orig" coincides with the View_center or if "dir" is
1480     parallel to r formed by projecting "orig" on the sphere. In
1481     either case, do a single test against the cell containing the
1482     intersection of "dir" and the sphere
1483     */
1484     /* orig will be updated-so preserve original value */
1485     if(!smMesh)
1486     return;
1487     point_on_sphere(b,orig,SM_VIEW_CENTER(smMesh));
1488     d = -DOT(b,dir);
1489     if(EQUAL_VEC3(orig,SM_VIEW_CENTER(smMesh)) || EQUAL(fabs(d),1.0))
1490     {
1491 gwlarson 3.8 qt = smPointLocateCell(smMesh,dir);
1492 gwlarson 3.5 /* Test triangles in the set for intersection with Ray:returns
1493     first found
1494     */
1495 gwlarson 3.8 #ifdef DEBUG
1496     if(QT_IS_EMPTY(qt))
1497     {
1498     eputs("smFindSamp(): point not found");
1499     return;
1500     }
1501     #endif
1502 gwlarson 3.5 ts = qtqueryset(qt);
1503     s_id = intersect_tri_set(ts,orig,dir,p);
1504     #ifdef DEBUG_TEST_DRIVER
1505     VCOPY(Pick_point[0],p);
1506     #endif
1507     }
1508     else
1509     {
1510 gwlarson 3.8 OBJECT t_set[QT_MAXCSET + 1];
1511     RT_ARGS rt;
1512    
1513 gwlarson 3.5 /* Test each of the root triangles against point id */
1514     QT_CLEAR_SET(t_set);
1515     VSUB(o,orig,SM_VIEW_CENTER(smMesh));
1516     ST_CLEAR_FLAGS(SM_LOCATOR(smMesh));
1517 gwlarson 3.8 rt.t_id = -1;
1518     rt.os = t_set;
1519     VCOPY(rt.orig,orig);
1520     VCOPY(rt.dir,dir);
1521     stTrace_ray(SM_LOCATOR(smMesh),o,dir,ray_trace_check_set,&rt);
1522     s_id = rt.t_id;
1523 gwlarson 3.5 }
1524     return(s_id);
1525     }
1526    
1527    
1528    
1529    
1530    
1531    
1532    
1533    
1534    
1535 gwlarson 3.1
1536