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.4 by gwlarson, Fri Sep 11 11:52:27 1998 UTC vs.
Revision 3.11 by gwlarson, Fri Mar 5 16:32:50 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 "sm_list.h"
17 + #include "sm_flag.h"
18   #include "sm_geom.h"
19 + #include "sm_qtree.h"
20   #include "sm_stree.h"
21  
14 /* Define 4 vertices on the sphere to create a tetrahedralization on
15   the sphere: triangles are as follows:
16        (2,1,0),(3,2,0), (1,3,0), (2,3,1)
17 */
18
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 < FVECT stDefault_base[4] = { {SQRT3_INV, SQRT3_INV, SQRT3_INV},
27 <                          {-SQRT3_INV, -SQRT3_INV, SQRT3_INV},
28 <                          {-SQRT3_INV, SQRT3_INV, -SQRT3_INV},
29 <                          {SQRT3_INV, -SQRT3_INV, -SQRT3_INV}};
30 < int stTri_verts[4][3] = { {2,1,0},{3,2,0},{1,3,0},{2,3,1}};
31 < int stTri_nbrs[4][3] =  { {2,1,3},{0,2,3},{1,0,3},{2,0,1}};
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 < stNth_base_verts(st,i,v1,v2,v3)
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;
32 int i;
33 FVECT v1,v2,v3;
63   {
64 <  VCOPY(v1,ST_NTH_BASE(st,stTri_verts[i][0]));
65 <  VCOPY(v2,ST_NTH_BASE(st,stTri_verts[i][1]));
66 <  VCOPY(v3,ST_NTH_BASE(st,stTri_verts[i][2]));
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 4 quadtrees rooted at st */
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;
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 <  /* stree always has 4 children corresponding to the base tris
98 <  */
99 <  for (i = 0; i < 4; i++)
100 <    qtFree(ST_NTH_ROOT(st, i));
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 <  QT_CLEAR_CHILDREN(ST_ROOT(st));
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 < STREE
129 < *stInit(st,center,base)
130 < STREE *st;
131 < FVECT  center,base[4];
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(base)
138 <    ST_SET_BASE(st,base);
137 >  if(STR_NTH_INDEX(root,0)==-1)
138 >    d0 = -v[0];
139    else
140 <    ST_SET_BASE(st,stDefault_base);
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 <  ST_SET_CENTER(st,center);
151 <  stClear(st);
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 <  return(st);
156 >  BARY_INT(d0,b[0])
157 >  BARY_INT(d1,b[1])
158 >  BARY_INT(d2,b[2])
159   }
160  
161  
74 /* "base" defines 4 vertices on the sphere to create a tetrahedralization on
75     the sphere: triangles are as follows:(0,1,2),(0,2,3), (0,3,1), (1,3,2)
76     if base is null: does default.
162  
163 < */
164 < STREE
165 < *stAlloc(st)
166 < STREE *st;
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(!st)
175 <    st = (STREE *)malloc(sizeof(STREE));
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 <  ST_ROOT(st) = qtAlloc();
211 <    
212 <  QT_CLEAR_CHILDREN(ST_ROOT(st));
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 <  return(st);
225 < }
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 < /* Find location of sample point in the DAG and return lowest level
237 <   containing triangle. "type" indicates whether the point was found
238 <   to be in interior to the triangle: GT_FACE, on one of its
239 <   edges GT_EDGE or coinciding with one of its vertices GT_VERTEX.
240 <   "which" specifies which vertex (0,1,2) or edge (0=v0v1, 1 = v1v2, 2 = v21)
241 < */
242 < int
243 < stPoint_locate(st,npt)
244 <    STREE *st;
105 <    FVECT npt;
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,d,j,id;
247 <    QUADTREE *rootptr,*qtptr;
109 <    FVECT v1,v2,v3;
110 <    OBJECT os[QT_MAXSET+1],*optr;
111 <    FVECT p0,p1,p2;
246 >  int i;
247 >  double d0,d1,d2;
248  
249 <    /* Test each of the root triangles against point id */
250 <    for(i=0; i < 4; i++)
251 <     {
252 <         rootptr = ST_NTH_ROOT_PTR(st,i);
253 <         stNth_base_verts(st,i,v1,v2,v3);
254 <         /* Return tri that p falls in */
255 <         qtptr = qtRoot_point_locate(rootptr,v1,v2,v3,npt,NULL,NULL,NULL);
256 <         if(!qtptr || QT_IS_EMPTY(*qtptr))
257 <            continue;
258 <         /* Get the set */
259 <         optr = qtqueryset(*qtptr);
260 <         for (j = QT_SET_CNT(optr),optr = QT_SET_PTR(optr);j > 0; j--)
261 <         {
262 <             /* Find the first triangle that pt falls */
263 <             id = QT_SET_NEXT_ELEM(optr);
264 <             qtTri_from_id(id,NULL,NULL,NULL,p0,p1,p2,NULL,NULL,NULL);
129 <             d = point_in_stri(p0,p1,p2,npt);  
130 <             if(d)
131 <               return(id);
132 <         }
133 <     }
134 <    return(EMPTY);
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_cell(st,p,t0,t1,t2)
270 > stPoint_locate(st,p)
271      STREE *st;
272      FVECT p;
141    FVECT t0,t1,t2;
273   {
274 <    int i,d;
275 <    QUADTREE *rootptr,*qtptr;
276 <    FVECT v0,v1,v2;
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 <    /* Test each of the root triangles against point id */
283 <    for(i=0; i < 4; i++)
284 <     {
285 <         rootptr = ST_NTH_ROOT_PTR(st,i);
286 <         stNth_base_verts(st,i,v0,v1,v2);
287 <         /* Return quadtree tri that p falls in */
288 <         qtptr = qtRoot_point_locate(rootptr,v0,v1,v2,p,t0,t1,t2);
289 <         if(qtptr)
290 <            return(qtptr);
291 <     }    /* Point not found */
292 <    return(NULL);
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 <
298 < int
299 < stAdd_tri(st,id,v0,v1,v2)
300 < STREE *st;
301 < int id;
302 < FVECT v0,v1,v2;
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 <  
305 <  int i,found;
306 <  QUADTREE *rootptr;
307 <  FVECT t0,t1,t2;
304 >  unsigned int cells,cross;
305 >  unsigned int vcell[3];
306 >  double t0,t1;
307 >  int i,inext;
308  
309 <  found = 0;
310 <  for(i=0; i < 4; i++)
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 <    rootptr = ST_NTH_ROOT_PTR(st,i);
322 <    stNth_base_verts(st,i,t0,t1,t2);
323 <    found |= qtRoot_add_tri(rootptr,t0,t1,t2,v0,v1,v2,id,0);
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(found);
374 >  return(cells);
375   }
376  
377 < int
378 < stApply_to_tri_cells(st,v0,v1,v2,func,arg)
379 < STREE *st;
380 < FVECT v0,v1,v2;
381 < int (*func)();
382 < int *arg;
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 <  int i,found;
387 <  QUADTREE *rootptr;
388 <  FVECT t0,t1,t2;
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 <  found = 0;
394 <  func(ST_ROOT_PTR(st),arg);
395 <  QT_SET_FLAG(ST_ROOT(st));
396 <  for(i=0; i < 4; i++)
397 <  {
199 <    rootptr = ST_NTH_ROOT_PTR(st,i);
200 <    stNth_base_verts(st,i,t0,t1,t2);
201 <    found |= qtApply_to_tri_cells(rootptr,v0,v1,v2,t0,t1,t2,func,arg);
202 <  }
203 <  return(found);
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 < int
452 < stRemove_tri(st,id,v0,v1,v2)
451 >
452 > stVisit_poly(st,verts,l,root,func,n)
453   STREE *st;
454 < int id;
455 < FVECT v0,v1,v2;
454 > FVECT *verts;
455 > LIST *l;
456 > unsigned int root;
457 > FUNC func;
458 > int n;
459   {
460 <  
461 <  int i,found;
217 <  QUADTREE *rootptr;
218 <  FVECT t0,t1,t2;
460 >  int id0,id1,id2;
461 >  FVECT tri[3];
462  
463 <  found = 0;
464 <  for(i=0; i < 4; i++)
463 >  id0 = pop_list(&l);
464 >  id1 = pop_list(&l);
465 >  while(l)
466    {
467 <    rootptr = ST_NTH_ROOT_PTR(st,i);
468 <    stNth_base_verts(st,i,t0,t1,t2);
469 <   found |= qtRemove_tri(rootptr,id,v0,v1,v2,t0,t1,t2);
467 >    id2 = pop_list(&l);
468 >    VCOPY(tri[0],verts[id0]);
469 >    VCOPY(tri[1],verts[id1]);
470 >    VCOPY(tri[2],verts[id2]);
471 >    stRoot_visit_tri(st,root,tri,func,n);
472 >    id1 = id2;
473    }
227  return(found);
474   }
475  
476 < int
477 < stVisit_tri_edges(st,t0,t1,t2,func,arg1,arg2)
478 <   STREE *st;
479 <   FVECT t0,t1,t2;
480 <   int (*func)();
481 <   int *arg1,arg2;
476 > stVisit_clip(st,i,verts,vcnt,l,cell,func,n)
477 >     STREE *st;
478 >     int i;
479 >     FVECT *verts;
480 >     int *vcnt;
481 >     LIST *l;
482 >     unsigned int cell;
483 >     FUNC func;
484 >     int n;
485   {
237    int id,i,w;
238    QUADTREE *rootptr;
239    FVECT q0,q1,q2,n,v[3],sdir[3],dir[3],tv,d;
240    double pd,t;
486  
487 <    VCOPY(v[0],t0); VCOPY(v[1],t1); VCOPY(v[2],t2);
488 <    VSUB(dir[0],t1,t0); VSUB(dir[1],t2,t1);VSUB(dir[2],t0,t2);
489 <    VCOPY(sdir[0],dir[0]);VCOPY(sdir[1],dir[1]);VCOPY(sdir[2],dir[2]);
490 <    w = 0;
491 <    for(i=0; i < 4; i++)
247 <     {
248 < #ifdef TEST_DRIVER
249 < Pick_cnt = 0;
250 < #endif
251 <         rootptr = ST_NTH_ROOT_PTR(st,i);
252 <         stNth_base_verts(st,i,q0,q1,q2);
253 <         /* Return quadtree tri that p falls in */
254 <         if(!point_in_stri(q0,q1,q2,v[w]))  
255 <           continue;
256 <         id = qtRoot_visit_tri_edges(rootptr,q0,q1,q2,v,dir,&w,func,arg1,arg2);
257 <         if(id == INVALID)
258 <         {
259 < #ifdef DEBUG
260 <           eputs("stVisit_tri_edges(): Unable to trace edges\n");
261 < #endif
262 <           return(INVALID);
263 <         }
264 <         if(id == QT_DONE)
265 <            return(*arg1);
266 <        
267 <         /* Crossed over to next cell: id = nbr */
268 <         while(1)
269 <         {
270 <             /* test if ray crosses plane between this quadtree triangle and
271 <                its neighbor- if it does then find intersection point with
272 <                ray and plane- this is the new origin
273 <                */
274 <           if(id==0)
275 <             VCROSS(n,q1,q2);
276 <           else
277 <             if(id==1)
278 <               VCROSS(n,q2,q0);
279 <           else
280 <             VCROSS(n,q0,q1);
487 >  LIST *labove,*lbelow,*endb,*enda;
488 >  int last = -1;
489 >  int id,last_id;
490 >  int first,first_id;
491 >  unsigned int cellb;
492  
493 <           if(w==0)
494 <             VCOPY(tv,t0);
495 <           else
496 <             if(w==1)
497 <               VCOPY(tv,t1);
498 <           else
499 <             VCOPY(tv,t2);
500 <           if(!intersect_ray_plane(tv,sdir[w],n,0.0,&t,v[w]))
501 <             return(INVALID);
493 >  labove = lbelow = NULL;
494 >  enda = endb = NULL;
495 >  while(l)
496 >  {
497 >    id = pop_list(&l);
498 >    if(ZERO(verts[id][i]))
499 >    {
500 >      if(last==-1)
501 >      {/* add below and above */
502 >        first = 2;
503 >        first_id= id;
504 >      }
505 >      lbelow=add_data(lbelow,id,&endb);
506 >      labove=add_data(labove,id,&enda);
507 >      last_id = id;
508 >      last = 2;
509 >      continue;
510 >    }
511 >    if(verts[id][i] < 0)
512 >    {
513 >      if(last != 1)
514 >      {
515 >        lbelow=add_data(lbelow,id,&endb);
516 >        if(last==-1)
517 >        {
518 >          first = 0;
519 >          first_id = id;
520 >        }
521 >        last_id = id;
522 >        last = 0;
523 >        continue;
524 >      }
525 >      /* intersect_edges */
526 >      intersect_edge_coord_plane(verts[last_id],verts[id],i,verts[*vcnt]);
527 >      /*newpoint goes to above and below*/
528 >      lbelow=add_data(lbelow,*vcnt,&endb);
529 >      lbelow=add_data(lbelow,id,&endb);
530 >      labove=add_data(labove,*vcnt,&enda);
531 >      last = 0;
532 >      last_id = id;
533 >      (*vcnt)++;
534 >    }
535 >    else
536 >    {
537 >      if(last != 0)
538 >      {
539 >        labove=add_data(labove,id,&enda);
540 >        if(last==-1)
541 >        {
542 >          first = 1;
543 >          first_id = id;
544 >        }
545 >        last_id = id;
546 >        last = 1;
547 >        continue;
548 >      }
549 >      /* intersect_edges */
550 >      /*newpoint goes to above and below*/
551 >      intersect_edge_coord_plane(verts[last_id],verts[id],i,verts[*vcnt]);
552 >      lbelow=add_data(lbelow,*vcnt,&endb);
553 >      labove=add_data(labove,*vcnt,&enda);
554 >      labove=add_data(labove,id,&enda);
555 >      last_id = id;
556 >      (*vcnt)++;
557 >      last = 1;
558 >    }
559 >  }
560 >  if(first != 2 && first != last)
561 >  {
562 >    intersect_edge_coord_plane(verts[id],verts[first_id],i,verts[*vcnt]);
563 >    /*newpoint goes to above and below*/
564 >    lbelow=add_data(lbelow,*vcnt,&endb);
565 >    labove=add_data(labove,*vcnt,&enda);
566 >    (*vcnt)++;
567  
568 <           VSUM(v[w],v[w],sdir[w],10.0*FTINY);
568 >  }
569 >  if(i==2)
570 >  {
571 >    if(lbelow)
572 >    {
573 >      if(LIST_NEXT(lbelow) && LIST_NEXT(LIST_NEXT(lbelow)))
574 >      {
575 >        cellb = cell | (1 << i);
576 >        stVisit_poly(st,verts,lbelow,cellb,func,n);
577 >      }
578 >      else
579 >        free_list(lbelow);
580 >    }
581 >    if(labove)
582 >     {
583 >      if(LIST_NEXT(labove) && LIST_NEXT(LIST_NEXT(labove)))
584 >        stVisit_poly(st,verts,labove,cell,func,n);
585 >      else
586 >        free_list(labove);
587 >     }
588 >  }
589 >  else
590 >  {
591 >    if(lbelow)
592 >    {
593 >      if(LIST_NEXT(lbelow) && LIST_NEXT(LIST_NEXT(lbelow)))
594 >        {
595 >          cellb = cell | (1 << i);
596 >          stVisit_clip(st,i+1,verts,vcnt,lbelow,cellb,func,n);
597 >        }
598 >      else
599 >        free_list(lbelow);
600 >    }
601 >    if(labove)
602 >     {
603 >       if(LIST_NEXT(labove) && LIST_NEXT(LIST_NEXT(labove)))
604 >         stVisit_clip(st,i+1,verts,vcnt,labove,cell,func,n);
605 >       else
606 >         free_list(labove);
607 >     }
608 >  }
609  
294           t = (1.0-t-10.0*FTINY);
295           if(t <= 0.0)
296           {
297             t = FTINY;
298 #if 0
299             eputs("stVisit_tri_edges(): edge end on plane\n");
300 #endif
301           }
302           dir[w][0] = sdir[w][0] * t;
303           dir[w][1] = sdir[w][1] * t;
304           dir[w][2] = sdir[w][2] * t;
305           i = stTri_nbrs[i][id];
306           rootptr = ST_NTH_ROOT_PTR(st,i);
307           stNth_base_verts(st,i,q0,q1,q2);
308           id=qtRoot_visit_tri_edges(rootptr,q0,q1,q2,v,dir,&w,func,arg1,arg2);
309           if(id == QT_DONE)
310             return(*arg1);
311           if(id == INVALID)
312             {
313 #if 0
314             eputs("stVisit_tri_edges(): point not found\n");
315 #endif  
316             return(INVALID);
317             }
318          
319         }
320     }    /* Point not found */
321    return(INVALID);
610   }
611  
612 <
325 < int
326 < stVisit_tri_edges2(st,t0,t1,t2,func,arg1,arg2)
612 > stVisit(st,tri,func,n)
613     STREE *st;
614 <   FVECT t0,t1,t2;
615 <   int (*func)();
616 <   int *arg1,arg2;
614 >   FVECT tri[3];
615 >   FUNC func;
616 >   int n;
617   {
618 <    int id,i,w;
619 <    QUADTREE *rootptr;
334 <    FVECT q0,q1,q2,v[3],i_pt;
618 >    int r0,r1,r2;
619 >    LIST *l;
620  
621 <    VCOPY(v[0],t0); VCOPY(v[1],t1); VCOPY(v[2],t2);
622 <    w = -1;
623 <    for(i=0; i < 4; i++)
624 <     {
625 < #ifdef TEST_DRIVER
626 < Pick_cnt = 0;
627 < #endif
628 <         rootptr = ST_NTH_ROOT_PTR(st,i);
629 <         stNth_base_verts(st,i,q0,q1,q2);
630 <         /* Return quadtree tri that p falls in */
631 <         if(!point_in_stri(q0,q1,q2,v[0]))  
632 <           continue;
633 <         id = qtRoot_visit_tri_edges2(rootptr,q0,q1,q2,v,i_pt,&w,
349 <                                      func,arg1,arg2);
350 <         if(id == INVALID)
351 <         {
352 < #ifdef DEBUG
353 <           eputs("stVisit_tri_edges(): Unable to trace edges\n");
354 < #endif
355 <           return(INVALID);
356 <         }
357 <         if(id == QT_DONE)
358 <            return(*arg1);
621 >    r0 = stLocate_root(tri[0]);
622 >    r1 = stLocate_root(tri[1]);
623 >    r2 = stLocate_root(tri[2]);
624 >    if(r0 == r1 && r1==r2)
625 >      stRoot_visit_tri(st,r0,tri,func,n);
626 >    else
627 >      {
628 >        FVECT verts[ST_CLIP_VERTS];
629 >        int cnt;
630 >
631 >        VCOPY(verts[0],tri[0]);
632 >        VCOPY(verts[1],tri[1]);
633 >        VCOPY(verts[2],tri[2]);
634          
635 <         /* Crossed over to next cell: id = nbr */
636 <         while(1)
637 <         {
638 <             /* test if ray crosses plane between this quadtree triangle and
639 <                its neighbor- if it does then find intersection point with
640 <                ray and plane- this is the new origin
366 <                */
367 <           i = stTri_nbrs[i][id];
368 <           rootptr = ST_NTH_ROOT_PTR(st,i);
369 <           stNth_base_verts(st,i,q0,q1,q2);
370 <           id=qtRoot_visit_tri_edges2(rootptr,q0,q1,q2,v,i_pt,&w,
371 <                                      func,arg1,arg2);
372 <           if(id == QT_DONE)
373 <             return(*arg1);
374 <           if(id == INVALID)
375 <             {
376 < #ifdef DEBUG
377 <             eputs("stVisit_tri_edges(): point not found\n");
378 < #endif  
379 <             return(INVALID);
380 <             }
381 <          
382 <         }
383 <     }    /* Point not found */
384 <    return(INVALID);
635 >        l = add_data(NULL,0,NULL);
636 >        l = add_data(l,1,NULL);
637 >        l = add_data(l,2,NULL);
638 >        cnt = 3;
639 >        stVisit_clip(st,0,verts,&cnt,l,0,func,n);
640 >      }
641   }
642  
387 int
388 stTrace_edge(st,orig,dir,max_t,func,arg1,arg2)
389   STREE *st;
390   FVECT orig,dir;
391   double max_t;
392   int (*func)();
393   int *arg1,arg2;
394 {
395    int id,i;
396    QUADTREE *rootptr;
397    FVECT q0,q1,q2,o,n,d;
398    double pd,t;
643  
644 < #if DEBUG
401 <    if(max_t > 1.0 || max_t < 0.0)
402 <    {
403 <      eputs("stTrace_edge(): max_t must be in [0,1]:adjusting\n");
404 <      max_t = 1.0;
405 <    }
406 < #endif
644 > /* New Insertion code!!! */
645  
408    VCOPY(o,orig);
409    for(i=0; i < 4; i++)
410     {
411 #ifdef TEST_DRIVER
412 Pick_cnt = 0;
413 #endif
414         rootptr = ST_NTH_ROOT_PTR(st,i);
415         stNth_base_verts(st,i,q0,q1,q2);
416         /* Return quadtree tri that p falls in */
417         id= qtRoot_trace_edge(rootptr,q0,q1,q2,o,dir,max_t,func,arg1,arg2);
418         if(id == INVALID)
419           continue;
420         if(id == QT_DONE)
421            return(*arg1);
422        
423         /* Crossed over to next cell: id = nbr */
424         while(1)
425         {
426             /* test if ray crosses plane between this quadtree triangle and
427                its neighbor- if it does then find intersection point with
428                ray and plane- this is the new origin
429                */
430           if(id==0)
431             VCROSS(n,q1,q2);
432           else
433             if(id==1)
434               VCROSS(n,q2,q0);
435           else
436             VCROSS(n,q0,q1);
646  
647 <           /* Ray does not cross into next cell: done and tri not found*/
439 <           if(!intersect_ray_plane(orig,dir,n,0.0,&t,o))
440 <             return(INVALID);
647 > BCOORD qtRoot[3][3] = { {MAXBCOORD2,0,0},{0,MAXBCOORD2,0},{0,0,MAXBCOORD2}};
648  
442           VSUM(o,o,dir,10*FTINY);
649  
650 <           d[0] = dir[0]*(1-t-10*FTINY);
651 <           d[1] = dir[1]*(1-t-10*FTINY);
652 <           d[2] = dir[2]*(1-t-10*FTINY);
653 <           i = stTri_nbrs[i][id];
654 <           rootptr = ST_NTH_ROOT_PTR(st,i);
655 <           stNth_base_verts(st,i,q0,q1,q2);
656 <           id = qtRoot_trace_edge(rootptr,q0,q1,q2,o,d,max_t,func,arg1,arg2);
657 <           if(id == QT_DONE)
658 <             return(*arg1);
659 <           if(id == INVALID)
660 <             {
661 < #if 0
662 <             eputs("stTrace_edges(): point not found\n");
663 < #endif  
664 <             return(INVALID);
459 <             }
460 <          
461 <         }
462 <     }    /* Point not found */
463 <    return(INVALID);
650 > convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02)
651 > int root;
652 > FVECT tri[3];
653 > BCOORD b0[3],b1[3],b2[3];
654 > BCOORD db10[3],db21[3],db02[3];
655 > {
656 >  /* Project the vertex into the qtree plane */
657 >  vert_to_qt_frame(root,tri[0],b0);
658 >  vert_to_qt_frame(root,tri[1],b1);
659 >  vert_to_qt_frame(root,tri[2],b2);
660 >
661 >  /* calculate triangle edge differences in new frame */
662 >  db10[0] = b1[0] - b0[0]; db10[1] = b1[1] - b0[1]; db10[2] = b1[2] - b0[2];
663 >  db21[0] = b2[0] - b1[0]; db21[1] = b2[1] - b1[1]; db21[2] = b2[2] - b1[2];
664 >  db02[0] = b0[0] - b2[0]; db02[1] = b0[1] - b2[1]; db02[2] = b0[2] - b2[2];
665   }
666  
667  
668 <
669 < int
469 < stTrace_ray(st,orig,dir,func,arg1,arg2)
668 > QUADTREE
669 > stRoot_insert_tri(st,root,tri,f)
670     STREE *st;
671 <   FVECT orig,dir;
672 <   int (*func)();
673 <   int *arg1,arg2;
671 >   int root;
672 >   FVECT tri[3];
673 >   FUNC f;
674   {
675 <    int id,i;
676 <    QUADTREE *rootptr;
677 <    FVECT q0,q1,q2,o,n;
678 <    double pd,t;
675 >  BCOORD b0[3],b1[3],b2[3];
676 >  BCOORD db10[3],db21[3],db02[3];
677 >  unsigned int s0,s1,s2,sq0,sq1,sq2;
678 >  QUADTREE qt;
679  
680 <    VCOPY(o,orig);
681 <    for(i=0; i < 4; i++)
482 <     {
483 < #ifdef TEST_DRIVER
484 < Pick_cnt = 0;
485 < #endif
486 <         rootptr = ST_NTH_ROOT_PTR(st,i);
487 <         stNth_base_verts(st,i,q0,q1,q2);
488 <         /* Return quadtree tri that p falls in */
489 <         id= qtRoot_trace_ray(rootptr,q0,q1,q2,o,dir,func,arg1,arg2);
490 <         if(id == INVALID)
491 <           continue;
492 <         if(id == QT_DONE)
493 <            return(*arg1);
494 <        
495 <         /* Crossed over to next cell: id = nbr */
496 <         while(1)
497 <         {
498 <             /* test if ray crosses plane between this quadtree triangle and
499 <                its neighbor- if it does then find intersection point with
500 <                ray and plane- this is the new origin
501 <                */
502 <           if(id==0)
503 <             VCROSS(n,q1,q2);
504 <           else
505 <             if(id==1)
506 <               VCROSS(n,q2,q0);
507 <           else
508 <             VCROSS(n,q0,q1);
680 >  /* Map the triangle vertices into the canonical barycentric frame */
681 >  convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02);
682  
683 <           /* Ray does not cross into next cell: done and tri not found*/
684 <           if(!intersect_ray_plane(orig,dir,n,0.0,NULL,o))
685 <             return(INVALID);
683 >  /* Calculate initial sidedness info */
684 >  SIDES_GTR(b0,b1,b2,s0,s1,s2,qtRoot[1][0],qtRoot[0][1],qtRoot[0][2]);
685 >  SIDES_GTR(b0,b1,b2,sq0,sq1,sq2,qtRoot[0][0],qtRoot[1][1],qtRoot[2][2]);
686  
687 <           VSUM(o,o,dir,10*FTINY);
688 <           i = stTri_nbrs[i][id];
689 <           rootptr = ST_NTH_ROOT_PTR(st,i);
690 <           stNth_base_verts(st,i,q0,q1,q2);
691 <           id = qtRoot_trace_ray(rootptr,q0,q1,q2,o,dir,func,arg1,arg2);
692 <           if(id == QT_DONE)
520 <             return(*arg1);
521 <           if(id == INVALID)
522 <             return(INVALID);
523 <          
524 <         }
525 <     }    /* Point not found */
526 <    return(INVALID);
687 >  qt = ST_ROOT_QT(st,root);
688 >  /* Visit cells that triangle intersects */
689 >  qt = qtInsert_tri(root,qt,qtRoot[0],qtRoot[1],qtRoot[2],
690 >       b0,b1,b2,db10,db21,db02,MAXBCOORD2 >> 1,s0,s1,s2, sq0,sq1,sq2,f,0);
691 >
692 >  return(qt);
693   }
694  
695 + stRoot_visit_tri(st,root,tri,f,n)
696 +   STREE *st;
697 +   int root;
698 +   FVECT tri[3];
699 +   FUNC f;
700 +   int n;
701 + {
702 +  BCOORD b0[3],b1[3],b2[3];
703 +  BCOORD db10[3],db21[3],db02[3];
704 +  unsigned int s0,s1,s2,sq0,sq1,sq2;
705 +  QUADTREE qt;
706  
707 +  /* Map the triangle vertices into the canonical barycentric frame */
708 +  convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02);
709  
710 < stVisit_tri_interior(st,t0,t1,t2,func,arg1,arg2)
710 >  /* Calculate initial sidedness info */
711 >  SIDES_GTR(b0,b1,b2,s0,s1,s2,qtRoot[1][0],qtRoot[0][1],qtRoot[0][2]);
712 >  SIDES_GTR(b0,b1,b2,sq0,sq1,sq2,qtRoot[0][0],qtRoot[1][1],qtRoot[2][2]);
713 >
714 >  qt = ST_ROOT_QT(st,root);
715 >  QT_SET_FLAG(ST_QT(st,root));
716 >  /* Visit cells that triangle intersects */
717 >  qtVisit_tri(root,qt,qtRoot[0],qtRoot[1],qtRoot[2],
718 >       b0,b1,b2,db10,db21,db02,MAXBCOORD2 >> 1,s0,s1,s2, sq0,sq1,sq2,f,n);
719 >
720 > }
721 >
722 > stInsert_tri(st,tri,f)
723     STREE *st;
724 <   FVECT t0,t1,t2;
725 <   int (*func)();
535 <   int *arg1,arg2;
724 >   FVECT tri[3];
725 >   FUNC f;
726   {
727 <  int i;
728 <  QUADTREE *rootptr;
729 <  FVECT q0,q1,q2;
727 >  unsigned int cells,which;
728 >  int root;
729 >  
730  
731 <  for(i=0; i < 4; i++)
731 >  /* calculate entry/exit points of edges through the cells */
732 >  cells = stTri_cells(st,tri);
733 >
734 >  /* For each cell that quadtree intersects: Map the triangle vertices into
735 >     the canonical barycentric frame of (1,0,0), (0,1,0),(0,0,1). Insert
736 >     by first doing a trivial reject on the interior nodes, and then a
737 >     tri/tri intersection at the leaf nodes.
738 >  */
739 >  for(root=0,which=1; root < ST_NUM_ROOT_NODES; root++,which <<= 1)
740    {
741 <    rootptr = ST_NTH_ROOT_PTR(st,i);
742 <    stNth_base_verts(st,i,q0,q1,q2);
743 <    qtVisit_tri_interior(rootptr,q0,q1,q2,t0,t1,t2,0,func,arg1,arg2);
741 >    /* For each of the quadtree roots: check if marked as intersecting tri*/
742 >    if(cells & which)
743 >      /* Visit tri cells */
744 >      ST_ROOT_QT(st,root) = stRoot_insert_tri(st,root,tri,f);
745    }
746   }
747  
748  
550 int
551 stApply_to_tri(st,t0,t1,t2,func,arg1,arg2)
552   STREE *st;
553   FVECT t0,t1,t2;
554   int (*func)();
555   int *arg1,arg2;
556 {
557    int f;
558    FVECT dir;
559    
560  /* First add all of the leaf cells lying on the triangle perimeter:
561     mark all cells seen on the way
562   */
563    qtClearAllFlags();          /* clear all quadtree branch flags */
564    f = 0;
565    VSUB(dir,t1,t0);
566    stTrace_edge(st,t0,dir,1.0,func,arg1,arg2);
567    VSUB(dir,t2,t1);
568    stTrace_edge(st,t1,dir,1.0,func,arg1,arg2);
569    VSUB(dir,t0,t2);
570    stTrace_edge(st,t2,dir,1.0,func,arg1,arg2);
571    /* Now visit interior */
572    stVisit_tri_interior(st,t0,t1,t2,func,arg1,arg2);
573 }
749  
750  
751  
752  
753  
579 int
580 stUpdate_tri(st,t_id,t0,t1,t2,edge_func,interior_func)
581   STREE *st;
582   int t_id;
583   FVECT t0,t1,t2;
584   int (*edge_func)(),(*interior_func)();
585 {
586    int f;
587    FVECT dir;
588    
589  /* First add all of the leaf cells lying on the triangle perimeter:
590     mark all cells seen on the way
591   */
592    ST_CLEAR_FLAGS(st);
593    f = 0;
594    /* Visit cells along edges of the tri */
595
596    stVisit_tri_edges2(st,t0,t1,t2,edge_func,&f,t_id);
597
598    /* Now visit interior */
599    if(QT_FLAG_FILL_TRI(f) || QT_FLAG_UPDATE(f))
600       stVisit_tri_interior(st,t0,t1,t2,interior_func,&f,t_id);
601 }
754  
755  
756  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines