--- ray/src/hd/sm_qtree.c	1999/03/05 16:32:49	3.12
+++ ray/src/hd/sm_qtree.c	2003/06/30 14:59:12	3.17
@@ -1,9 +1,6 @@
-/* Copyright (c) 1998 Silicon Graphics, Inc. */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ SGI";
+static const char	RCSid[] = "$Id: sm_qtree.c,v 3.17 2003/06/30 14:59:12 schorsch Exp $";
 #endif
-
 /*
 *  sm_qtree.c: adapted from octree.c from radiance code
 */
@@ -13,25 +10,30 @@ static char SCCSid[] = "$SunId$ SGI";
  *     7/28/85
  */
 
+#include <string.h>
+
 #include "standard.h"
 #include "sm_flag.h"
 #include "sm_geom.h"
+#include "object.h"
 #include "sm_qtree.h"
 
-QUADTREE  *quad_block[QT_MAX_BLK];		/* our quadtree */
+QUADTREE  *quad_block[QT_MAX_BLK];	  /* our quadtree */
 static QUADTREE  quad_free_list = EMPTY;  /* freed octree nodes */
-static QUADTREE  treetop = 0;		/* next free node */
-int4 *quad_flag= NULL;
+static QUADTREE  treetop = 0;		  /* next free node */
+int32 *quad_flag= NULL;                    /* quadtree flags */
 
-#ifdef TEST_DRIVER
-extern FVECT Pick_v0[500],Pick_v1[500],Pick_v2[500];
-extern int Pick_cnt,Pick_tri,Pick_samp;
-extern FVECT Pick_point[500];
-extern int Pick_q[500];
 
+qtremovelast(qt,id)
+QUADTREE qt;
+OBJECT id;
+{
+#ifdef DEBUG
+  if(qtqueryset(qt)[(*(qtqueryset(qt)))] != id)
+    error(CONSISTENCY,"not removed\n");
 #endif
-int Incnt=0;
-
+  ((*(qtqueryset(qt)))--);
+}
 QUADTREE
 qtAlloc()			/* allocate a quadtree */
 {
@@ -52,7 +54,7 @@ qtAlloc()			/* allocate a quadtree */
 	   error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
 
 	/* Realloc the per/node flags */
-	quad_flag = (int4 *)realloc((char *)quad_flag,
+	quad_flag = (int32 *)realloc((void *)quad_flag,
 			(QT_BLOCK(freet)+1)*((QT_BLOCK_SIZE+7)>>3));
 	if (quad_flag == NULL)
 	   error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
@@ -68,7 +70,8 @@ qtClearAllFlags()		/* clear all quadtree branch flags 
     return;
   
   /* Clear the node flags*/
-  bzero((char *)quad_flag, (QT_BLOCK(treetop-4)+1)*((QT_BLOCK_SIZE+7)>>3));
+  memset((char *)quad_flag, '\0',
+		  (QT_BLOCK(treetop-4)+1)*((QT_BLOCK_SIZE+7)>>3));
   /* Clear set flags */
   qtclearsetflags();
 }
@@ -99,16 +102,173 @@ qtDone()			/* free EVERYTHING */
 	{
 	    if (quad_block[i] == NULL)
 		break;
-	    free((char *)quad_block[i]);
+	    free((void *)quad_block[i]);
 	    quad_block[i] = NULL;
 	}
 	/* Free the flags */
-	if (i) free((char *)quad_flag);
+	if (i) free((void *)quad_flag);
 	quad_flag = NULL;
 	quad_free_list = EMPTY;
 	treetop = 0;
 }
 
+/*
+ * bary_parent(coord,i)  : Set coord to be value of pt at one level up in
+ *                         quadtree, given it was the ith child
+ * BCOORD coord[3];      :barycentric coordinates of point, also will hold
+ *                        result as side effect
+ * int i;                :designates which child coord was
+ */
+bary_parent(coord,i)
+BCOORD coord[3];
+int i;
+{
+  switch(i) {
+  case 0:
+    /* update bary for child */
+    coord[0] = (coord[0] >> 1) + MAXBCOORD4;
+    coord[1] >>= 1;
+    coord[2] >>= 1;
+    break;
+  case 1:
+    coord[0] >>= 1;
+    coord[1]  = (coord[1] >> 1) + MAXBCOORD4;
+    coord[2] >>= 1;
+    break;
+    
+  case 2:
+    coord[0] >>= 1;
+    coord[1] >>= 1;
+    coord[2] = (coord[2] >> 1) + MAXBCOORD4;
+    break;
+    
+  case 3:
+    coord[0] = MAXBCOORD4 - (coord[0] >> 1);
+    coord[1] = MAXBCOORD4 - (coord[1] >> 1);
+    coord[2] = MAXBCOORD4 - (coord[2] >> 1);
+    break;
+#ifdef DEBUG
+  default:
+    eputs("bary_parent():Invalid child\n");
+    break;
+#endif
+  }
+}
+
+/*
+ * bary_from_child(coord,child,next): Given that coord was the (child) child
+ *                      of the quadtree node, calculate what the (next) child
+ *                      is at the same level.
+ * BCOORD coord[3];   :At current level (child)th child of node,will also hold
+ *                     result as side effect
+ * int child,next;    :Which child coord was (child), and which one
+ *                     is now desired(next)
+ */
+bary_from_child(coord,child,next)
+BCOORD coord[3];
+int child,next;
+{
+#ifdef DEBUG
+  if(child <0 || child > 3)
+  {
+    eputs("bary_from_child():Invalid child\n");
+    return;
+  }
+  if(next <0 || next > 3)
+  {
+    eputs("bary_from_child():Invalid next\n");
+    return;
+  }
+#endif
+  if(next == child)
+    return;
+
+  switch(child){
+  case 0:
+      coord[0] = 0;
+      coord[1] = MAXBCOORD2 - coord[1];
+      coord[2] = MAXBCOORD2 - coord[2];
+      break;
+  case 1:
+      coord[0] = MAXBCOORD2 - coord[0];
+      coord[1] = 0;
+      coord[2] = MAXBCOORD2 - coord[2];
+      break;
+  case 2:
+      coord[0] = MAXBCOORD2 - coord[0];
+      coord[1] = MAXBCOORD2 - coord[1];
+      coord[2] = 0;
+    break;
+  case 3:
+    switch(next){
+    case 0:
+      coord[0] = 0;
+      coord[1] = MAXBCOORD2 - coord[1];
+      coord[2] = MAXBCOORD2 - coord[2];
+      break;
+    case 1:
+      coord[0] = MAXBCOORD2 - coord[0];
+      coord[1] = 0;
+      coord[2] = MAXBCOORD2 - coord[2];
+      break;
+    case 2:
+      coord[0] = MAXBCOORD2 - coord[0];
+      coord[1] = MAXBCOORD2 - coord[1];
+      coord[2] = 0;
+      break;
+    }
+    break;
+  }
+}
+
+/*
+ * int
+ * bary_child(coord):  Return which child coord is of node at current level
+ *                     As side effect recalculate coord at next level
+ * BCOORD coord[3];   :At current level coordinates of point, will also set
+ *                     to coords at next level as side effect
+ */
+int
+bary_child(coord)
+BCOORD coord[3];
+{
+  int i;
+
+  if(coord[0] > MAXBCOORD4)
+  { 
+      /* update bary for child */
+      coord[0] = (coord[0]<< 1) - MAXBCOORD2;
+      coord[1] <<= 1;
+      coord[2] <<= 1;
+      return(0);
+  }
+  else
+    if(coord[1] > MAXBCOORD4)
+    {
+      coord[0] <<= 1;
+      coord[1] = (coord[1] << 1) - MAXBCOORD2;
+      coord[2] <<= 1;
+      return(1);
+    }
+    else
+      if(coord[2] > MAXBCOORD4)
+      {
+	coord[0] <<= 1;
+	coord[1] <<= 1;
+	coord[2] = (coord[2] << 1) - MAXBCOORD2;
+	return(2);
+      }
+      else
+	 {
+	   coord[0] = MAXBCOORD2 - (coord[0] << 1);
+	   coord[1] = MAXBCOORD2 - (coord[1] << 1);
+	   coord[2] = MAXBCOORD2 - (coord[2] << 1);
+	   return(3);
+	 }
+}
+
+
+
 QUADTREE
 qtLocate(qt,bcoord)
 QUADTREE qt;
@@ -247,13 +407,6 @@ qtTrace_ray(qt,b,db0,db1,db2,nextptr,sign,sfactor,func
    }
 }    
 
-
-
-
-
-
-
-
 #define TEST_INT(tl,th,d,q0,q1,h,l) \
                   tl=d>q0;th=d<q1; if(tl&&th) goto Lfunc_call; \
 		 if(tl) if(l) goto Lfunc_call; else h=TRUE; \
@@ -281,18 +434,24 @@ int n;
 
   /* First check if can trivial accept: if vertex in cell */
   if(s0 & s1 & s2)
+  {
     goto Lfunc_call;
-
+  }
   /* Assumption: Could not trivial reject*/
   /* IF any coords lie on edge- must be in if couldnt reject */
   a = q1[0];b= q0[1]; c= q0[2];
   if(t0[0]== a || t0[1] == b || t0[2] == c)
+  {
     goto Lfunc_call;
+  }
   if(t1[0]== a || t1[1] == b || t1[2] == c)
+  {
     goto Lfunc_call;
+  }
   if(t2[0]== a || t2[1] == b || t2[2] == c)
+  {
     goto Lfunc_call;
-  
+  }
   /* Test for edge crossings */
   /* Test t0t1,t1t2,t2t0 against a */
   /* Calculate edge crossings */
@@ -411,26 +570,33 @@ int n;
   char testl,testh,test_t0t1,test_t1t2,test_t2t0;
   unsigned int ls0,ls1,ls2;
   
+  
   /* May have gotten passed trivial reject if one/two vertices are ON */
   a = q1[0];b= q0[1]; c= q0[2];
   SIDES_LESS(t0,t1,t2,ls0,ls1,ls2,a,b,c);
   
   /* First check if can trivial accept: if vertex in cell */
   if(ls0 & ls1 & ls2)
+  {
     goto Lfunc_call;
-
+  }
   if(ls0==0 || ls1 == 0 || ls2 ==0)
     return(qt);
   /* Assumption: Could not trivial reject*/
   /* IF any coords lie on edge- must be in if couldnt reject */
 
   if(t0[0]== a || t0[1] == b || t0[2] == c)
+  {
     goto Lfunc_call;
+  }
   if(t1[0]== a || t1[1] == b || t1[2] == c)
+  {
     goto Lfunc_call;
+  }
   if(t2[0]== a || t2[1] == b || t2[2] == c)
+  {
     goto Lfunc_call;
-
+  }
   /* Test for edge crossings */
   /* Test t0t1 against a,b,c */
   /* First test if t0t1 can be trivially rejected */
@@ -454,7 +620,7 @@ int n;
       db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
       TEST_INT(testl,testh,db,q1[1],b,bl,bh)
   }
-  test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) || 
+ test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) || 
        ((sq0 & 5)==0) ||((sq1 & 5)==0)|| ((sq2 & 5)==0));
   if(test_t2t0 && (e0 & 4))
   {
@@ -519,7 +685,6 @@ int n;
 
   return(qt);
 Lfunc_call:
-
   qt = f.func(f.argptr,root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
 	      s0,s1,s2,sq0,sq1,sq2,1,f,n);
   return(qt);
@@ -650,6 +815,7 @@ BCOORD dt10[3],dt21[3],dt02[3];
 BCOORD scale;
 unsigned int s0,s1,s2,sq0,sq1,sq2;
 FUNC f;
+int n;
 {
   BCOORD a,b,c;
   BCOORD va[3],vb[3],vc[3];
@@ -740,7 +906,7 @@ FUNC f;
 /*************************************************************************/
 /* Visit code for applying functions: NOTE THIS IS THE SAME AS INSERT CODE:
   except sets flags: wanted insert to be as efficient as possible: so 
-  broke into two sets of routines
+  broke into two sets of routines. Also traverses only to n levels.
 */
 
 qtVisit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
@@ -1205,6 +1371,195 @@ Lfunc_call:
   if(!QT_IS_EMPTY(qt))
     QT_LEAF_SET_FLAG(qt);
 }
+
+
+QUADTREE
+qtInsert_point(root,qt,parent,q0,q1,q2,bc,scale,f,n,doneptr)
+int root;
+QUADTREE qt,parent;
+BCOORD q0[3],q1[3],q2[3],bc[3],scale; 
+FUNC f;
+int n,*doneptr;
+{
+  BCOORD a,b,c;
+  BCOORD va[3],vb[3],vc[3];
+
+  if(QT_IS_TREE(qt))
+  {  
+    a = q1[0] + scale;
+    b = q0[1] + scale;
+    c = q0[2] + scale;
+    if(bc[0] > a)
+    {
+      vc[0] = a; vc[1] = b; vc[2] = q0[2];
+      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
+      QT_NTH_CHILD(qt,0) = qtInsert_point(root,QT_NTH_CHILD(qt,0),
+				  qt,q0,vc,vb,bc,scale>>1,f,n+1,doneptr);
+      if(!*doneptr)
+      {
+	f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
+      }	      
+      return(qt);
+    }
+    if(bc[1] > b)
+    {
+      vc[0] = a; vc[1] = b; vc[2] = q0[2];
+      va[0] = q1[0]; va[1] = b; va[2] = c;
+      QT_NTH_CHILD(qt,1) =qtInsert_point(root,QT_NTH_CHILD(qt,1),
+				 qt,vc,q1,va,bc,scale >>1,f,n+1,doneptr);
+      if(!*doneptr)
+      {
+	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
+      }	      
+      return(qt);
+    }
+    if(bc[2] > c)
+    {
+      va[0] = q1[0]; va[1] = b; va[2] = c;
+      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
+      QT_NTH_CHILD(qt,2) = qtInsert_point(root,QT_NTH_CHILD(qt,2),
+				  qt,vb,va,q2,bc,scale>>1,f,n+1,doneptr);
+      if(!*doneptr)
+      {
+	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
+      }	
+      return(qt);
+    }
+    else
+    {
+      va[0] = q1[0]; va[1] = b; va[2] = c;
+      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
+      vc[0] = a; vc[1] = b; vc[2] = q0[2];
+      QT_NTH_CHILD(qt,3) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,3),
+				      qt,va,vb,vc,bc,scale>>1,f,n+1,doneptr);
+      if(!*doneptr)
+      {
+	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
+      }	
+      return(qt);
+    }
+  }
+  else
+  {
+    qt = f.func(f.argptr,root,qt,parent,q0,q1,q2,bc,scale,0,f,n,doneptr);
+    return(qt);
+  }
+}
+
+
+QUADTREE
+qtInsert_point_rev(root,qt,parent,q0,q1,q2,bc,scale,f,n,doneptr)
+int root;
+QUADTREE qt,parent;
+BCOORD q0[3],q1[3],q2[3],bc[3]; 
+BCOORD scale;
+FUNC f;
+int n,*doneptr;
+{
+  BCOORD a,b,c;
+  BCOORD va[3],vb[3],vc[3];
+
+  if(QT_IS_TREE(qt))
+  {  
+    a = q1[0] - scale;
+    b = q0[1] - scale;
+    c = q0[2] - scale;
+    if(bc[0] < a)
+    {
+      vc[0] = a; vc[1] = b; vc[2] = q0[2];
+      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
+      QT_NTH_CHILD(qt,0) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,0),
+				  qt,q0,vc,vb,bc,scale>>1,f,n+1,doneptr);
+      if(!*doneptr)
+      {
+	f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
+      }	
+      return(qt);
+    }
+    if(bc[1] < b)
+    {
+      vc[0] = a; vc[1] = b; vc[2] = q0[2];
+      va[0] = q1[0]; va[1] = b; va[2] = c;
+      QT_NTH_CHILD(qt,1) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,1),
+				   qt,vc,q1,va,bc,scale>>1,f,n+1,doneptr);
+      if(!*doneptr)
+      {
+	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
+      }	
+      return(qt);
+    }
+    if(bc[2] < c)
+    {
+      va[0] = q1[0]; va[1] = b; va[2] = c;
+      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
+      QT_NTH_CHILD(qt,2) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,2),
+				qt,vb,va,q2,bc,scale>>1,f,n+1,doneptr);
+      if(!*doneptr)
+      {
+	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
+      }	
+      return(qt);
+    }
+    else
+    {
+      va[0] = q1[0]; va[1] = b; va[2] = c;
+      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
+      vc[0] = a; vc[1] = b; vc[2] = q0[2];
+      QT_NTH_CHILD(qt,3) = qtInsert_point(root,QT_NTH_CHILD(qt,3),
+				   qt,va,vb,vc,bc,scale>>1,f,n+1,doneptr);
+      if(!*doneptr)
+      {
+	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
+	if(!*doneptr)
+	  f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
+      }	
+      return(qt);
+    }
+  }
+  else
+  {
+    qt = f.func(f.argptr,root,qt,parent,q0,q1,q2,bc,scale,1,f,n,doneptr);
+    return(qt);
+  }
+}
+
+
+
+
+
+
+
+