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.13 by greg, Sat Feb 22 02:07:25 2003 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines