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

Comparing ray/src/hd/sm_stree.c (file contents):
Revision 3.3 by gwlarson, Tue Aug 25 11:03:27 1998 UTC vs.
Revision 3.12 by gwlarson, Thu Jun 10 15:22:24 1999 UTC

# Line 6 | Line 6 | static char SCCSid[] = "$SunId$ SGI";
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 "object.h"
17 <
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 < /* Define 4 vertices on the sphere to create a tetrahedralization on
24 <   the sphere: triangles are as follows:
25 <        (0,1,2),(0,2,3), (0,3,1), (1,3,2)
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  
22 FVECT stDefault_base[4] = { {SQRT3_INV, SQRT3_INV, SQRT3_INV},
23                          {-SQRT3_INV, -SQRT3_INV, SQRT3_INV},
24                          {-SQRT3_INV, SQRT3_INV, -SQRT3_INV},
25                          {SQRT3_INV, -SQRT3_INV, -SQRT3_INV}};
26 int stTri_verts[4][3] = { {2,1,0},
27                          {3,2,0},
28                          {1,3,0},
29                          {2,3,1}};
61  
62 < stNth_base_verts(st,i,v1,v2,v3)
62 > stInit(st)
63   STREE *st;
33 int i;
34 FVECT v1,v2,v3;
64   {
65 <  VCOPY(v1,ST_NTH_BASE(st,stTri_verts[i][0]));
66 <  VCOPY(v2,ST_NTH_BASE(st,stTri_verts[i][1]));
67 <  VCOPY(v3,ST_NTH_BASE(st,stTri_verts[i][2]));
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 < /* Frees the 4 quadtrees rooted at st */
42 < stClear(st)
77 > stFree(st)
78   STREE *st;
79   {
80 <  int i;
80 >  qtDone();
81 >  free(st);
82 > }
83  
84 <  /* stree always has 4 children corresponding to the base tris
85 <  */
86 <  for (i = 0; i < 4; i++)
87 <    qtFree(ST_NTH_ROOT(st, i));
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 <  QT_CLEAR_CHILDREN(ST_ROOT(st));
98 <
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 < STREE
107 < *stInit(st,center,base)
108 < STREE *st;
109 < FVECT  center,base[4];
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(base)
116 <    ST_SET_BASE(st,base);
115 >  if(STR_NTH_INDEX(root,0)==-1)
116 >    d0 = -v[0];
117    else
118 <    ST_SET_BASE(st,stDefault_base);
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 <  ST_SET_CENTER(st,center);
129 <  stClear(st);
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 <  return(st);
134 >  BARY_INT(d0,b[0])
135 >  BARY_INT(d1,b[1])
136 >  BARY_INT(d2,b[2])
137   }
138  
139  
75 /* "base" defines 4 vertices on the sphere to create a tetrahedralization on
76     the sphere: triangles are as follows:(0,1,2),(0,2,3), (0,3,1), (1,3,2)
77     if base is null: does default.
140  
141 < */
142 < STREE
143 < *stAlloc(st)
144 < STREE *st;
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(!st)
153 <    st = (STREE *)malloc(sizeof(STREE));
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 <  ST_ROOT(st) = qtAlloc();
189 <    
190 <  QT_CLEAR_CHILDREN(ST_ROOT(st));
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 <  return(st);
203 < }
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 < /* Find location of sample point in the DAG and return lowest level
215 <   containing triangle. "type" indicates whether the point was found
216 <   to be in interior to the triangle: GT_FACE, on one of its
217 <   edges GT_EDGE or coinciding with one of its vertices GT_VERTEX.
218 <   "which" specifies which vertex (0,1,2) or edge (0=v0v1, 1 = v1v2, 2 = v21)
219 < */
220 < int
221 < stPoint_locate(st,npt,type,which)
222 <    STREE *st;
106 <    FVECT npt;
107 <    char *type,*which;
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,d,j,id;
225 <    QUADTREE *rootptr,*qtptr;
111 <    FVECT v1,v2,v3;
112 <    OBJECT os[MAXSET+1],*optr;
113 <    char w;
114 <    FVECT p0,p1,p2;
224 >  int i;
225 >  double d0,d1,d2;
226  
227 <    /* Test each of the root triangles against point id */
228 <    for(i=0; i < 4; i++)
229 <     {
230 <         rootptr = ST_NTH_ROOT_PTR(st,i);
231 <         stNth_base_verts(st,i,v1,v2,v3);
232 <         /* Return tri that p falls in */
233 <         qtptr = qtRoot_point_locate(rootptr,v1,v2,v3,npt,NULL,NULL,NULL);
234 <         if(!qtptr)
235 <            continue;
236 <         /* Get the set */
237 <         qtgetset(os,*qtptr);
238 <         for (j = QT_SET_CNT(os),optr = QT_SET_PTR(os); j > 0; j--)
239 <         {
240 <             /* Find the first triangle that pt falls */
241 <             id = QT_SET_NEXT_ELEM(optr);
242 <             qtTri_verts_from_id(id,p0,p1,p2);
132 <             d = test_single_point_against_spherical_tri(p0,p1,p2,npt,&w);  
133 <             if(d)
134 <                {
135 <                    if(type)
136 <                       *type = d;
137 <                    if(which)
138 <                       *which = w;
139 <                    return(id);
140 <                }
141 <         }
142 <     }
143 <    if(which)
144 <      *which = 0;
145 <    if(type)
146 <      *type = 0;
147 <    return(EMPTY);
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_cell(st,p,t0,t1,t2)
248 > stPoint_locate(st,p)
249      STREE *st;
250      FVECT p;
154    FVECT t0,t1,t2;
251   {
252 <    int i,d;
253 <    QUADTREE *rootptr,*qtptr;
254 <    FVECT v0,v1,v2;
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 <    /* Test each of the root triangles against point id */
261 <    for(i=0; i < 4; i++)
262 <     {
263 <         rootptr = ST_NTH_ROOT_PTR(st,i);
264 <         stNth_base_verts(st,i,v0,v1,v2);
265 <         /* Return tri that p falls in */
266 <         qtptr = qtRoot_point_locate(rootptr,v0,v1,v2,p,t0,t1,t2);
267 <         /* NOTE: For now return only one triangle */
268 <         if(qtptr)
269 <            return(qtptr);
270 <     }    /* Point not found */
271 <    return(NULL);
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 < QUADTREE
288 < *stAdd_tri_from_pt(st,p,t_id)
289 <    STREE *st;
290 <    FVECT p;
291 <    int t_id;
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 <    int i,d;
365 <    QUADTREE *rootptr,*qtptr;
366 <    FVECT v0,v1,v2;
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 <    
372 <    /* Test each of the root triangles against point id */
373 <    for(i=0; i < 4; i++)
374 <     {
375 <         rootptr = ST_NTH_ROOT_PTR(st,i);
191 <         stNth_base_verts(st,i,v0,v1,v2);
192 <         /* Return tri that p falls in */
193 <         qtptr = qtRoot_add_tri_from_point(rootptr,v0,v1,v2,p,t_id);
194 <         /* NOTE: For now return only one triangle */
195 <         if(qtptr)
196 <            return(qtptr);
197 <     }    /* Point not found */
198 <    return(NULL);
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 < stAdd_tri(st,id,v0,v1,v2)
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 < int id;
433 < FVECT v0,v1,v2;
432 > FVECT *verts;
433 > LIST *l;
434 > unsigned int root;
435 > FUNC func;
436 > int n;
437   {
438 <  
439 <  int i,found;
209 <  QUADTREE *rootptr;
210 <  FVECT t0,t1,t2;
438 >  int id0,id1,id2;
439 >  FVECT tri[3];
440  
441 <  found = 0;
442 <  for(i=0; i < 4; i++)
441 >  id0 = pop_list(&l);
442 >  id1 = pop_list(&l);
443 >  while(l)
444    {
445 <    rootptr = ST_NTH_ROOT_PTR(st,i);
446 <    stNth_base_verts(st,i,t0,t1,t2);
447 <    found |= qtRoot_add_tri(rootptr,id,v0,v1,v2,t0,t1,t2,0);
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    }
219  return(found);
452   }
453  
454 < int
455 < stApply_to_tri_cells(st,v0,v1,v2,func,arg)
456 < STREE *st;
457 < FVECT v0,v1,v2;
458 < int (*func)();
459 < char *arg;
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   {
229  int i,found;
230  QUADTREE *rootptr;
231  FVECT t0,t1,t2;
464  
465 <  found = 0;
466 <  for(i=0; i < 4; i++)
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 <    rootptr = ST_NTH_ROOT_PTR(st,i);
476 <    stNth_base_verts(st,i,t0,t1,t2);
477 <    found |= qtApply_to_tri_cells(rootptr,v0,v1,v2,t0,t1,t2,func,arg);
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 <  return(found);
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 < int
622 < stRemove_tri(st,id,v0,v1,v2)
623 < STREE *st;
624 < int id;
625 < FVECT v0,v1,v2;
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 <  
635 <  int i,found;
636 <  QUADTREE *rootptr;
637 <  FVECT t0,t1,t2;
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 <  found = 0;
640 <  for(i=0; i < 4; i++)
641 <  {
642 <    rootptr = ST_NTH_ROOT_PTR(st,i);
261 <    stNth_base_verts(st,i,t0,t1,t2);
262 <   found |= qtRemove_tri(rootptr,id,v0,v1,v2,t0,t1,t2);
263 <  }
264 <  return(found);
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 < #if 0
666 < int
667 < stAdd_tri_opt(st,id,v0,v1,v2)
668 < STREE *st;
669 < int id;
670 < FVECT v0,v1,v2;
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    
279  int i,found;
280  QUADTREE *qtptr;
281  FVECT pt,t0,t1,t2;
708  
709 <  /* First add all of the leaf cells lying on the triangle perimeter:
710 <     mark all cells seen on the way
285 <   */
286 <  /* clear all of the flags */
287 <  qtClearAllFlags();            /* clear all quadtree branch flags */
709 >  /* calculate entry/exit points of edges through the cells */
710 >  cells = stTri_cells(st,tri);
711  
712 <  /* Now trace each triangle edge-marking cells visited, and adding tri to
713 <     the leafs
714 <   */
715 <  stAdd_tri_from_pt(st,v0,id,t0,t1,t2);
716 <  /* Find next cell that projection of ray intersects */
717 <  VCOPY(pt,v0);
295 <  /* NOTE: Check if in same cell */
296 <  while(traceEdge(pt,v1,t0,t1,t2,pt))
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 <      stAdd_tri_from_pt(st,pt,id,t0,t1,t2);
720 <      traceEdge(pt,v1,t0,t1,t2,pt);
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 <  while(traceEdge(pt,v2,t0,t1,t2,pt))
302 <  {
303 <      stAdd_tri_from_pt(st,pt,id,t0,t1,t2);
304 <      traceEdge(pt,v2,t0,t1,t2,pt);
305 <  }
306 <  while(traceEdge(pt,v0,t0,t1,t2,pt))
307 <  {
308 <      stAdd_tri_from_pt(st,pt,id,t0,t1,t2);
309 <      traceEdge(pt,v2,t0,t1,t2,pt);
310 <  }
724 > }
725  
726 <  /* NOTE: Optimization: if <= 2 cells added: dont need to fill */
727 <  /* Traverse: follow nodes with flag set or one vertex in triangle */
728 <  
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 < #endif
751 >
752 >
753 >
754 >
755 >
756 >
757 >
758 >
759 >

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines