ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/sm_stree.c
Revision: 3.12
Committed: Thu Jun 10 15:22:24 1999 UTC (24 years, 10 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.11: +34 -32 lines
Log Message:
Implemented sample quadtree in place of triangle quadtree
Made geometric predicates more robust
Added #define LORES which utilizes a single precision floating point
  sample array, the default is a double sample array
Added topology DEBUG commands (for DEBUG > 1)
Made code optimizations

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