ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/sm_stree.c
Revision: 3.9
Committed: Tue Jan 5 16:52:39 1999 UTC (25 years, 3 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.8: +0 -0 lines
Log Message:
fixed logic in smTest to handle nearby and coincident points, added base tris to rendering, made list of new triangles to speed up rendering

File Contents

# Content
1 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ SGI";
5 #endif
6
7 /*
8 * sm_stree.c
9 * An stree (spherical quadtree) is defined by an octahedron in
10 * canonical form,and a world center point. Each face of the
11 * octahedron is adaptively subdivided as a planar triangular quadtree.
12 * World space geometry is projected onto the quadtree faces from the
13 * sphere center.
14 */
15 #include "standard.h"
16 #include "sm_list.h"
17 #include "sm_flag.h"
18 #include "sm_geom.h"
19 #include "sm_qtree.h"
20 #include "sm_stree.h"
21
22 #ifdef TEST_DRIVER
23 extern FVECT Pick_point[500],Pick_v0[500],Pick_v1[500],Pick_v2[500];
24 extern int Pick_cnt;
25 #endif
26 /* octahedron coordinates */
27 FVECT stDefault_base[6] = { {1.,0.,0.},{0.,1.,0.}, {0.,0.,1.},
28 {-1.,0.,0.},{0.,-1.,0.},{0.,0.,-1.}};
29 /* octahedron triangle vertices */
30 int stBase_verts[8][3] = { {0,1,2},{3,1,2},{0,4,2},{3,4,2},
31 {0,1,5},{3,1,5},{0,4,5},{3,4,5}};
32 /* octahedron triangle nbrs ; nbr i is the face opposite vertex i*/
33 int stBase_nbrs[8][3] = { {1,2,4},{0,3,5},{3,0,6},{2,1,7},
34 {5,6,0},{4,7,1},{7,4,2},{6,5,3}};
35 int stRoot_indices[8][3] = {{1,1,1},{-1,1,1},{1,-1,1},{-1,-1,1},
36 {1,1,-1},{-1,1,-1},{1,-1,-1},{-1,-1,-1}};
37 /*
38 +z y -z y
39 | |
40 1 | 0 5 | 4
41 ______|______ x _______|______ x
42 3 | 2 7 | 6
43 | |
44
45 Nbrs
46 +z y -z y
47 /0|1\ /1|0\
48 5 / | \ 4 / | \
49 /(1)|(0)\ 1 /(5)|(4)\ 0
50 / | \ / | \
51 /2 1|0 2\ /2 0|1 2\
52 /------|------\x /------|------\x
53 \0 1|2 0/ \0 2|2 1/
54 \ | / \ | /
55 7\ (3)|(2) / 6 3 \ (7)|(6) / 2
56 \ | / \ | /
57 \ 2|1 / \ 1|0 /
58 */
59
60
61 stInit(st)
62 STREE *st;
63 {
64 int i,j;
65
66 ST_TOP_QT(st) = qtAlloc();
67 ST_BOTTOM_QT(st) = qtAlloc();
68 /* Clear the children */
69
70 QT_CLEAR_CHILDREN(ST_TOP_QT(st));
71 QT_CLEAR_CHILDREN(ST_BOTTOM_QT(st));
72 }
73
74 /* Frees the children of the 2 quadtrees rooted at st,
75 Does not free root nodes: just clears
76 */
77 stClear(st)
78 STREE *st;
79 {
80 qtDone();
81 stInit(st);
82 }
83
84 /* Allocates a stree structure and creates octahedron base */
85 STREE
86 *stAlloc(st)
87 STREE *st;
88 {
89 int i,m;
90 FVECT v0,v1,v2;
91 FVECT n;
92
93 if(!st)
94 if(!(st = (STREE *)malloc(sizeof(STREE))))
95 error(SYSTEM,"stAlloc(): Unable to allocate memory\n");
96
97 /* Allocate the top and bottom quadtree root nodes */
98 stInit(st);
99
100
101 /* will go ********************************************/
102 /* Set the octahedron base */
103 ST_SET_BASE(st,stDefault_base);
104
105 /* Calculate octahedron face and edge normals */
106 for(i=0; i < ST_NUM_ROOT_NODES; i++)
107 {
108 VCOPY(v0,ST_NTH_V(st,i,0));
109 VCOPY(v1,ST_NTH_V(st,i,1));
110 VCOPY(v2,ST_NTH_V(st,i,2));
111 tri_plane_equation(v0,v1,v2, &ST_NTH_PLANE(st,i),FALSE);
112 m = max_index(FP_N(ST_NTH_PLANE(st,i)),NULL);
113 FP_X(ST_NTH_PLANE(st,i)) = (m+1)%3;
114 FP_Y(ST_NTH_PLANE(st,i)) = (m+2)%3;
115 FP_Z(ST_NTH_PLANE(st,i)) = m;
116 VCROSS(ST_EDGE_NORM(st,i,0),v0,v1);
117 VCROSS(ST_EDGE_NORM(st,i,1),v1,v2);
118 VCROSS(ST_EDGE_NORM(st,i,2),v2,v0);
119 }
120
121 /*****************************************************************/
122 return(st);
123 }
124
125 #define BARY_INT(v,b) if((v)>2.0) (b) = MAXBCOORD;else \
126 if((v)<-2.0) (b)=-MAXBCOORD;else (b)=(BCOORD)((v)*MAXBCOORD2);
127
128 vert_to_qt_frame(root,v,b)
129 int root;
130 FVECT v;
131 BCOORD b[3];
132 {
133 int i;
134 double scale;
135 double d0,d1,d2;
136
137 if(STR_NTH_INDEX(root,0)==-1)
138 d0 = -v[0];
139 else
140 d0 = v[0];
141 if(STR_NTH_INDEX(root,1)==-1)
142 d1 = -v[1];
143 else
144 d1 = v[1];
145 if(STR_NTH_INDEX(root,2)==-1)
146 d2 = -v[2];
147 else
148 d2 = v[2];
149
150 /* Plane is now x+y+z = 1 - intersection of pt ray is qtv/den */
151 scale = 1.0/ (d0 + d1 + d2);
152 d0 *= scale;
153 d1 *= scale;
154 d2 *= scale;
155
156 BARY_INT(d0,b[0])
157 BARY_INT(d1,b[1])
158 BARY_INT(d2,b[2])
159 }
160
161
162
163
164 ray_to_qt_frame(root,v,dir,b,db)
165 int root;
166 FVECT v,dir;
167 BCOORD b[3],db[3];
168 {
169 int i;
170 double scale;
171 double d0,d1,d2;
172 double dir0,dir1,dir2;
173
174 if(STR_NTH_INDEX(root,0)==-1)
175 {
176 d0 = -v[0];
177 dir0 = -dir[0];
178 }
179 else
180 {
181 d0 = v[0];
182 dir0 = dir[0];
183 }
184 if(STR_NTH_INDEX(root,1)==-1)
185 {
186 d1 = -v[1];
187 dir1 = -dir[1];
188 }
189 else
190 {
191 d1 = v[1];
192 dir1 = dir[1];
193 }
194 if(STR_NTH_INDEX(root,2)==-1)
195 {
196 d2 = -v[2];
197 dir2 = -dir[2];
198 }
199 else
200 {
201 d2 = v[2];
202 dir2 = dir[2];
203 }
204 /* Plane is now x+y+z = 1 - intersection of pt ray is qtv/den */
205 scale = 1.0/ (d0 + d1 + d2);
206 d0 *= scale;
207 d1 *= scale;
208 d2 *= scale;
209
210 /* Calculate intersection point of orig+dir: This calculation is done
211 after the origin is projected into the plane in order to constrain
212 the projection( i.e. the size of the projection of the unit direction
213 vector translated to the origin depends on how close
214 the origin is to the view center
215 */
216 /* Must divide by at least root2 to insure that projection will fit
217 int [-2,2] bounds: assumed length is 1: therefore greatest projection
218 from endpoint of triangle is at 45 degrees or projected length of root2
219 */
220 dir0 = d0 + dir0*0.5;
221 dir1 = d1 + dir1*0.5;
222 dir2 = d2 + dir2*0.5;
223
224 scale = 1.0/ (dir0 + dir1 + dir2);
225 dir0 *= scale;
226 dir1 *= scale;
227 dir2 *= scale;
228
229 BARY_INT(d0,b[0])
230 BARY_INT(d1,b[1])
231 BARY_INT(d2,b[2])
232 BARY_INT(dir0,db[0])
233 BARY_INT(dir1,db[1])
234 BARY_INT(dir2,db[2])
235
236 db[0] -= b[0];
237 db[1] -= b[1];
238 db[2] -= b[2];
239 }
240
241 qt_frame_to_vert(root,b,v)
242 int root;
243 BCOORD b[3];
244 FVECT v;
245 {
246 int i;
247 double d0,d1,d2;
248
249 d0 = b[0]/(double)MAXBCOORD2;
250 d1 = b[1]/(double)MAXBCOORD2;
251 d2 = b[2]/(double)MAXBCOORD2;
252
253 if(STR_NTH_INDEX(root,0)==-1)
254 v[0] = -d0;
255 else
256 v[0] = d0;
257 if(STR_NTH_INDEX(root,1)==-1)
258 v[1] = -d1;
259 else
260 v[1] = d1;
261 if(STR_NTH_INDEX(root,2)==-1)
262 v[2] = -d2;
263 else
264 v[2] = d2;
265 }
266
267
268 /* Return quadtree leaf node containing point 'p'*/
269 QUADTREE
270 stPoint_locate(st,p)
271 STREE *st;
272 FVECT p;
273 {
274 QUADTREE qt;
275 BCOORD bcoordi[3];
276 int i;
277
278 /* Find root quadtree that contains p */
279 i = stLocate_root(p);
280 qt = ST_ROOT_QT(st,i);
281
282 /* Will return lowest level triangle containing point: It the
283 point is on an edge or vertex: will return first associated
284 triangle encountered in the child traversal- the others can
285 be derived using triangle adjacency information
286 */
287 if(QT_IS_TREE(qt))
288 {
289 vert_to_qt_frame(i,p,bcoordi);
290 i = bary_child(bcoordi);
291 return(qtLocate(QT_NTH_CHILD(qt,i),bcoordi));
292 }
293 else
294 return(qt);
295 }
296
297 static unsigned int nbr_b[8][3] ={{2,4,16},{1,8,32},{8,1,64},{4,2,128},
298 {32,64,1},{16,128,2},{128,16,4},{64,32,8}};
299 unsigned int
300 stTri_cells(st,v)
301 STREE *st;
302 FVECT v[3];
303 {
304 unsigned int cells,cross;
305 unsigned int vcell[3];
306 double t0,t1;
307 int i,inext;
308
309 /* First find base cells that tri vertices are in (0-7)*/
310 vcell[0] = stLocate_root(v[0]);
311 vcell[1] = stLocate_root(v[1]);
312 vcell[2] = stLocate_root(v[2]);
313
314 /* If all in same cell- return that bit only */
315 if(vcell[0] == vcell[1] && vcell[1] == vcell[2])
316 return( 1 << vcell[0]);
317
318 cells = 0;
319 for(i=0;i<3; i++)
320 {
321 if(i==2)
322 inext = 0;
323 else
324 inext = i+1;
325 /* Mark cell containing initial vertex */
326 cells |= 1 << vcell[i];
327
328 /* Take the exclusive or: will have bits set where edge crosses axis=0*/
329 cross = vcell[i] ^ vcell[inext];
330 /* If crosses 2 planes: then have 2 options for edge crossing-pick closest
331 otherwise just hits two*/
332 /* Neighbors are zyx */
333 switch(cross){
334 case 3: /* crosses x=0 and y=0 */
335 t0 = -v[i][0]/(v[inext][0]-v[i][0]);
336 t1 = -v[i][1]/(v[inext][1]-v[i][1]);
337 if(t0==t1)
338 break;
339 else if(t0 < t1)
340 cells |= nbr_b[vcell[i]][0];
341 else
342 cells |= nbr_b[vcell[i]][1];
343 break;
344 case 5: /* crosses x=0 and z=0 */
345 t0 = -v[i][0]/(v[inext][0]-v[i][0]);
346 t1 = -v[i][2]/(v[inext][2]-v[i][2]);
347 if(t0==t1)
348 break;
349 else if(t0 < t1)
350 cells |= nbr_b[vcell[i]][0];
351 else
352 cells |=nbr_b[vcell[i]][2];
353
354 break;
355 case 6:/* crosses z=0 and y=0 */
356 t0 = -v[i][2]/(v[inext][2]-v[i][2]);
357 t1 = -v[i][1]/(v[inext][1]-v[i][1]);
358 if(t0==t1)
359 break;
360 else if(t0 < t1)
361 {
362 cells |= nbr_b[vcell[i]][2];
363 }
364 else
365 {
366 cells |=nbr_b[vcell[i]][1];
367 }
368 break;
369 case 7:
370 error(CONSISTENCY," Insert:Edge shouldnt be able to span 3 cells");
371 break;
372 }
373 }
374 return(cells);
375 }
376
377
378 stRoot_trace_ray(qt,root,orig,dir,nextptr,func,f)
379 QUADTREE qt;
380 int root;
381 FVECT orig,dir;
382 int *nextptr;
383 FUNC func;
384 int *f;
385 {
386 double br[3];
387 BCOORD bi[3],dbi[3];
388
389 /* Project the origin onto the root node plane */
390 /* Find the intersection point of the origin */
391 ray_to_qt_frame(root,orig,dir,bi,dbi);
392
393 /* trace the ray starting with this node */
394 qtTrace_ray(qt,bi,dbi[0],dbi[1],dbi[2],nextptr,0,0,func,f);
395 if(!QT_FLAG_IS_DONE(*f))
396 qt_frame_to_vert(root,bi,orig);
397
398 }
399
400 /* Trace ray 'orig-dir' through stree and apply 'func(qtptr,f,argptr)' at each
401 node that it intersects
402 */
403 int
404 stTrace_ray(st,orig,dir,func)
405 STREE *st;
406 FVECT orig,dir;
407 FUNC func;
408 {
409 int next,last,i,f=0;
410 QUADTREE qt;
411 FVECT o,n,v;
412 double pd,t,d;
413
414 VCOPY(o,orig);
415 #ifdef TEST_DRIVER
416 Pick_cnt=0;
417 #endif;
418 /* Find the qt node that o falls in */
419 i = stLocate_root(o);
420 qt = ST_ROOT_QT(st,i);
421
422 stRoot_trace_ray(qt,i,o,dir,&next,func,&f);
423
424 if(QT_FLAG_IS_DONE(f))
425 return(TRUE);
426 /*
427 d = DOT(orig,dir)/sqrt(DOT(orig,orig));
428 VSUM(v,orig,dir,-d);
429 */
430 /* Crossed over to next cell: id = nbr */
431 while(1)
432 {
433 /* test if ray crosses plane between this quadtree triangle and
434 its neighbor- if it does then find intersection point with
435 ray and plane- this is the new origin
436 */
437 if(next == INVALID)
438 return(FALSE);
439 /*
440 if(DOT(o,v) < 0.0)
441 return(FALSE);
442 */
443 i = stBase_nbrs[i][next];
444 qt = ST_ROOT_QT(st,i);
445 stRoot_trace_ray(qt,i,o,dir,&next,func,&f);
446 if(QT_FLAG_IS_DONE(f))
447 return(TRUE);
448 }
449 }
450
451
452 stVisit_poly(st,verts,l,root,func)
453 STREE *st;
454 FVECT *verts;
455 LIST *l;
456 unsigned int root;
457 FUNC func;
458 {
459 int id0,id1,id2;
460 FVECT tri[3];
461
462 id0 = pop_list(&l);
463 id1 = pop_list(&l);
464 while(l)
465 {
466 id2 = pop_list(&l);
467 VCOPY(tri[0],verts[id0]);
468 VCOPY(tri[1],verts[id1]);
469 VCOPY(tri[2],verts[id2]);
470 stRoot_visit_tri(st,root,tri,func);
471 id1 = id2;
472 }
473 }
474
475 stVisit_clip(st,i,verts,vcnt,l,cell,func)
476 STREE *st;
477 int i;
478 FVECT *verts;
479 int *vcnt;
480 LIST *l;
481 unsigned int cell;
482 FUNC func;
483 {
484
485 LIST *labove,*lbelow,*endb,*enda;
486 int last = -1;
487 int id,last_id;
488 int first,first_id;
489 unsigned int cellb;
490
491 labove = lbelow = NULL;
492 enda = endb = NULL;
493 while(l)
494 {
495 id = pop_list(&l);
496 if(ZERO(verts[id][i]))
497 {
498 if(last==-1)
499 {/* add below and above */
500 first = 2;
501 first_id= id;
502 }
503 lbelow=add_data(lbelow,id,&endb);
504 labove=add_data(labove,id,&enda);
505 last_id = id;
506 last = 2;
507 continue;
508 }
509 if(verts[id][i] < 0)
510 {
511 if(last != 1)
512 {
513 lbelow=add_data(lbelow,id,&endb);
514 if(last==-1)
515 {
516 first = 0;
517 first_id = id;
518 }
519 last_id = id;
520 last = 0;
521 continue;
522 }
523 /* intersect_edges */
524 intersect_edge_coord_plane(verts[last_id],verts[id],i,verts[*vcnt]);
525 /*newpoint goes to above and below*/
526 lbelow=add_data(lbelow,*vcnt,&endb);
527 lbelow=add_data(lbelow,id,&endb);
528 labove=add_data(labove,*vcnt,&enda);
529 last = 0;
530 last_id = id;
531 (*vcnt)++;
532 }
533 else
534 {
535 if(last != 0)
536 {
537 labove=add_data(labove,id,&enda);
538 if(last==-1)
539 {
540 first = 1;
541 first_id = id;
542 }
543 last_id = id;
544 last = 1;
545 continue;
546 }
547 /* intersect_edges */
548 /*newpoint goes to above and below*/
549 intersect_edge_coord_plane(verts[last_id],verts[id],i,verts[*vcnt]);
550 lbelow=add_data(lbelow,*vcnt,&endb);
551 labove=add_data(labove,*vcnt,&enda);
552 labove=add_data(labove,id,&enda);
553 last_id = id;
554 (*vcnt)++;
555 last = 1;
556 }
557 }
558 if(first != 2 && first != last)
559 {
560 intersect_edge_coord_plane(verts[id],verts[first_id],i,verts[*vcnt]);
561 /*newpoint goes to above and below*/
562 lbelow=add_data(lbelow,*vcnt,&endb);
563 labove=add_data(labove,*vcnt,&enda);
564 (*vcnt)++;
565
566 }
567 if(i==2)
568 {
569 if(lbelow)
570 {
571 if(LIST_NEXT(lbelow) && LIST_NEXT(LIST_NEXT(lbelow)))
572 {
573 cellb = cell | (1 << i);
574 stVisit_poly(st,verts,lbelow,cellb,func);
575 }
576 else
577 free_list(lbelow);
578 }
579 if(labove)
580 {
581 if(LIST_NEXT(labove) && LIST_NEXT(LIST_NEXT(labove)))
582 stVisit_poly(st,verts,labove,cell,func);
583 else
584 free_list(labove);
585 }
586 }
587 else
588 {
589 if(lbelow)
590 {
591 if(LIST_NEXT(lbelow) && LIST_NEXT(LIST_NEXT(lbelow)))
592 {
593 cellb = cell | (1 << i);
594 stVisit_clip(st,i+1,verts,vcnt,lbelow,cellb,func);
595 }
596 else
597 free_list(lbelow);
598 }
599 if(labove)
600 {
601 if(LIST_NEXT(labove) && LIST_NEXT(LIST_NEXT(labove)))
602 stVisit_clip(st,i+1,verts,vcnt,labove,cell,func);
603 else
604 free_list(labove);
605 }
606 }
607
608 }
609
610 stVisit(st,tri,func)
611 STREE *st;
612 FVECT tri[3];
613 FUNC func;
614 {
615 int r0,r1,r2;
616 LIST *l;
617
618 r0 = stLocate_root(tri[0]);
619 r1 = stLocate_root(tri[1]);
620 r2 = stLocate_root(tri[2]);
621 if(r0 == r1 && r1==r2)
622 stRoot_visit_tri(st,r0,tri,func);
623 else
624 {
625 FVECT verts[ST_CLIP_VERTS];
626 int cnt;
627
628 VCOPY(verts[0],tri[0]);
629 VCOPY(verts[1],tri[1]);
630 VCOPY(verts[2],tri[2]);
631
632 l = add_data(NULL,0,NULL);
633 l = add_data(l,1,NULL);
634 l = add_data(l,2,NULL);
635 cnt = 3;
636 stVisit_clip(st,0,verts,&cnt,l,0,func);
637 }
638 }
639
640
641 /* New Insertion code!!! */
642
643
644 BCOORD qtRoot[3][3] = { {MAXBCOORD2,0,0},{0,MAXBCOORD2,0},{0,0,MAXBCOORD2}};
645
646
647 convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02)
648 int root;
649 FVECT tri[3];
650 BCOORD b0[3],b1[3],b2[3];
651 BCOORD db10[3],db21[3],db02[3];
652 {
653 /* Project the vertex into the qtree plane */
654 vert_to_qt_frame(root,tri[0],b0);
655 vert_to_qt_frame(root,tri[1],b1);
656 vert_to_qt_frame(root,tri[2],b2);
657
658 /* calculate triangle edge differences in new frame */
659 db10[0] = b1[0] - b0[0]; db10[1] = b1[1] - b0[1]; db10[2] = b1[2] - b0[2];
660 db21[0] = b2[0] - b1[0]; db21[1] = b2[1] - b1[1]; db21[2] = b2[2] - b1[2];
661 db02[0] = b0[0] - b2[0]; db02[1] = b0[1] - b2[1]; db02[2] = b0[2] - b2[2];
662 }
663
664
665 QUADTREE
666 stRoot_insert_tri(st,root,tri,f)
667 STREE *st;
668 int root;
669 FVECT tri[3];
670 FUNC f;
671 {
672 BCOORD b0[3],b1[3],b2[3];
673 BCOORD db10[3],db21[3],db02[3];
674 unsigned int s0,s1,s2,sq0,sq1,sq2;
675 QUADTREE qt;
676
677 /* Map the triangle vertices into the canonical barycentric frame */
678 convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02);
679
680 /* Calculate initial sidedness info */
681 SIDES_GTR(b0,b1,b2,s0,s1,s2,qtRoot[1][0],qtRoot[0][1],qtRoot[0][2]);
682 SIDES_GTR(b0,b1,b2,sq0,sq1,sq2,qtRoot[0][0],qtRoot[1][1],qtRoot[2][2]);
683
684 qt = ST_ROOT_QT(st,root);
685 /* Visit cells that triangle intersects */
686 qt = qtInsert_tri(root,qt,qtRoot[0],qtRoot[1],qtRoot[2],
687 b0,b1,b2,db10,db21,db02,MAXBCOORD2 >> 1,s0,s1,s2, sq0,sq1,sq2,f,0);
688
689 return(qt);
690 }
691
692 stRoot_visit_tri(st,root,tri,f)
693 STREE *st;
694 int root;
695 FVECT tri[3];
696 FUNC f;
697 {
698 BCOORD b0[3],b1[3],b2[3];
699 BCOORD db10[3],db21[3],db02[3];
700 unsigned int s0,s1,s2,sq0,sq1,sq2;
701 QUADTREE qt;
702
703 /* Map the triangle vertices into the canonical barycentric frame */
704 convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02);
705
706 /* Calculate initial sidedness info */
707 SIDES_GTR(b0,b1,b2,s0,s1,s2,qtRoot[1][0],qtRoot[0][1],qtRoot[0][2]);
708 SIDES_GTR(b0,b1,b2,sq0,sq1,sq2,qtRoot[0][0],qtRoot[1][1],qtRoot[2][2]);
709
710 qt = ST_ROOT_QT(st,root);
711 QT_SET_FLAG(ST_QT(st,root));
712 /* Visit cells that triangle intersects */
713 qtVisit_tri(root,qt,qtRoot[0],qtRoot[1],qtRoot[2],
714 b0,b1,b2,db10,db21,db02,MAXBCOORD2 >> 1,s0,s1,s2, sq0,sq1,sq2,f);
715
716 }
717
718 stInsert_tri(st,tri,f)
719 STREE *st;
720 FVECT tri[3];
721 FUNC f;
722 {
723 unsigned int cells,which;
724 int root;
725
726
727 /* calculate entry/exit points of edges through the cells */
728 cells = stTri_cells(st,tri);
729
730 /* For each cell that quadtree intersects: Map the triangle vertices into
731 the canonical barycentric frame of (1,0,0), (0,1,0),(0,0,1). Insert
732 by first doing a trivial reject on the interior nodes, and then a
733 tri/tri intersection at the leaf nodes.
734 */
735 for(root=0,which=1; root < ST_NUM_ROOT_NODES; root++,which <<= 1)
736 {
737 /* For each of the quadtree roots: check if marked as intersecting tri*/
738 if(cells & which)
739 /* Visit tri cells */
740 ST_ROOT_QT(st,root) = stRoot_insert_tri(st,root,tri,f);
741 }
742 }
743
744
745
746
747
748
749
750
751
752
753