RAYTRACE.C

Go to the documentation of this file.

#define MODULE_TRACE

// This is the raytracing module - it uses an octree spatial
// decomposition to speed things up - see the book!

#include <setjmp.h>
#include "render.h"

#define UC unsigned short
#define US unsigned long

typedef struct traceNODE {
  long n,            
       mi[3],ma[3];  
  UC *vxi;           
  UC *oid;           
  US *fid;           
  struct traceNODE *octree[8]; 
  long id; struct traceNODE *last; long level; 
} node;

//#define THRESHOLDF 45     // nominal max number of faces per voxel
//#define THRESHOLDR 6      // maximum number of recursive levels   
#define THRESHOLDF 25     // nominal max number of faces per voxel - use when not subdividing 
#define THRESHOLDR 8      // maximum number of recursive levels  

extern short R_terminator(void);

static jmp_buf x_buf;
static long threshold_f,threshold_z,threshold_r,
            nfaces,kcount,maxfaces,minfaces,
            trace_depth_max=7, // number of reflections etc
            level_depth_max;
static node *ROOT;
static long *ObjectXlist=NULL;    // extra number of vertices required for this object due to voxel splitting
static long voxel_xmin,voxel_ymin,voxel_zmin,
            voxel_xmax,voxel_ymax,voxel_zmax;

static node *which_voxel(vector p, node *nc);
static node *track_ray(node *nc, vector p, vector d);
static void destroyBST(node *np);
static void countBST(node *np);
static void printBST(node *np);
static int MakeVoxels(node *np, int level);
static void eliminateDuplicatesInBST(node *n);
static BOOL PartitionVoxel(node *, node *, node *, long, long index);
static BOOL SplitVoxel(node *, node *, node *, long, long index);

static double longest_ray;

static void IntersectPlane(long p, long index, ivector vi1, ivector vi2, ivector vo){
 long id1[3]={1,2,0},id2[3]={2,0,1};
 double mu;
 mu=(double)(vi2[index]-vi1[index]);
 if(fabs(mu) < 1.0)mu=0.0;
 else mu=((double)p-vi1[index])/mu;
 vo[index]=p;
 vo[id1[index]]=vi1[id1[index]]+(long)(mu*((double)(vi2[id1[index]]-vi1[id1[index]])));
 vo[id2[index]]=vi1[id2[index]]+(long)(mu*((double)(vi2[id2[index]]-vi1[id2[index]])));
}

static double c_angle(ivector p1, ivector p2, ivector p3){ // angle at p2
 vector d1,d2;
 VECSUB(p1,p2,d1)
 VECSUB(p3,p2,d2)
 normalize(d1); normalize(d2);
 return DOT(d1,d2);
}

static void SplitFaceAtVoxelPlane(long o, long i, long c, // object face and type of intersection
                                  tvertex *v, tface *f,     // Base pointers 
                                  tvertex *v0, tvertex *v1, tvertex *v2, tface *fi,
                                  long p, long index,
                                  UC *oid1, US *fid1, UC *oid2, US *fid2,
                                  long *n1, long *n2){
 long vp,fp,fa,fb,i1,i2,i3,va,vb,i1i,i2i,i3i;
 tvertex *vl,*vr,*vrr;
 double uu,vv,ww,uu1,vv1,ww1,uui1,vvi1,wwi1,uui2,vvi2,wwi2,uui3,vvi3,wwi3;
 vector norm,normm,norm1,norm2,norm3;
 vp=Object[o].NoTraceVertices;
 fp=Object[o].NoTraceFaces;
 va=vp; vb=vp+1;
 fa=fp; fb=fp+1;
 //if(debug != NULL)fprintf(debug,"Splitting object %ld (v=%ld f=%ld) face %ld new vertices %ld %ld new faces %ld %ld\n",o,vp,fp,i,va,vb,fa,fb); 
 if(c == 3 || c == 4 || c == 5 || c == 13 || c == 14 || c == 15){  // add one vertex and one face
   if      (c ==  3 || c == 13){vl = v2; vr=v1; i1=0; i2=1; i3=2;}
   else if (c ==  4 || c == 14){vl = v0; vr=v2; i1=1; i2=2; i3=0;}
   else if (c ==  5 || c == 15){vl = v1; vr=v0; i1=2; i2=0; i3=1;}
   //if(debug != NULL)fprintf(debug,"Vertex in split plane !!\n");
   i1i=f[i].V[i1];  i2i=f[i].V[i2]; i3i=f[i].V[i3]; // actual vertex IDs
   IntersectPlane(p,index,vl->ip,vr->ip, v[va].ip);
   f[i].V[0]=i1i; f[i].V[1]=i2i; f[i].V[2]=va;
   f[fa].V[0]=i1i; f[fa].V[1]=va; f[fa].V[2]=i3i;
   f[fa].o=f[i].o; // new face points at same original face
   f[fa].f=f[i].f; // as the face that is being subdivided
   if(c == 3 || c == 4 || c == 5){
     oid1[*n1] = o; fid1[*n1] = fa; *n1 = (*n1) + 1;
     oid2[*n2] = o; fid2[*n2] =  i; *n2 = (*n2) + 1;
   }
   else{
     oid1[*n1] = o; fid1[*n1] =  i; *n1 = (*n1) + 1;
     oid2[*n2] = o; fid2[*n2] = fa; *n2 = (*n2) + 1;
   }
   vp += 1; fp += 1;
 }
 else {  // add two vertices and two faces 
   if     (c == 6 || c == 10){vl=v0; vr=v1; vrr=v2; i1=0; i2=1; i3=2;}
   else if(c == 7 || c == 11){vl=v1; vr=v2; vrr=v0; i1=1; i2=2; i3=0;}
   else if(c == 8 || c == 12){vl=v2; vr=v0; vrr=v1; i1=2; i2=0; i3=1;}
   i1i=f[i].V[i1];  i2i=f[i].V[i2]; i3i=f[i].V[i3]; // actual vertex IDs
   IntersectPlane(p,index,vl->ip,vr->ip, v[va].ip);
   IntersectPlane(p,index,vl->ip,vrr->ip, v[vb].ip);
   f[i].V[0]=i1i; f[i].V[1]=va; f[i].V[2]=vb;
   if(c_angle(v[i3i].ip,v[vb].ip,v[va].ip) >= c_angle(v[vb].ip,v[va].ip,v[i2i].ip)){
     f[fa].V[0]=i3i; f[fa].V[1]=vb; f[fa].V[2]=va;
     f[fb].V[0]=i2i; f[fb].V[1]=i3i; f[fb].V[2]=va;
   }
   else{
     f[fa].V[0]=i3i; f[fa].V[1]=vb; f[fa].V[2]=i2i;
     f[fb].V[0]=i2i; f[fb].V[1]=vb; f[fb].V[2]=va;
   }
   f[fa].o=f[i].o; // new face points at same original face
   f[fa].f=f[i].f; // as the face that is being subdivided
   f[fb].o=f[i].o; // new face points at same original face
   f[fb].f=f[i].f; // as the face that is being subdivided
   if(c == 6 || c == 7 || c == 8){
     oid1[*n1] = o; fid1[*n1] =  i; *n1 = (*n1) + 1;
     oid2[*n2] = o; fid2[*n2] = fa; *n2 = (*n2) + 1;
     oid2[*n2] = o; fid2[*n2] = fb; *n2 = (*n2) + 1;
   }
   else{
     oid1[*n1] = o; fid1[*n1] = fb; *n1 = (*n1) + 1;
     oid1[*n1] = o; fid1[*n1] = fa; *n1 = (*n1) + 1;
     oid2[*n2] = o; fid2[*n2] =  i; *n2 = (*n2) + 1;
   }
   //fprintf(debug,"X o=%ld i=%ld fa=%ld fb=%ld  n1=%ld n2=%ld\n",o,i,fa,fb,*n1,*n2);
   vp += 2; fp+=2;
 }
 Object[o].NoTraceVertices = vp;
 Object[o].NoTraceFaces = fp;
}

#define IN_PLANE(point,plane) (point == plane)
#define IN_LESS_VOXEL(point,plane) (point < plane)
#define IN_GREATER_VOXEL(point,plane) (point > plane)
#define IN_LESS_EQUAL_VOXEL(point,plane) (point <= plane)
#define IN_GREATER_EQUAL_VOXEL(point,plane) (point >= plane)

static long FaceIntersectVoxelPlane(tvertex *v0, tvertex *v1, tvertex *v2, long p, long index, UC *vxi){
 long p0,p1,p2;
 // determine the type of split needed for this face and return number of extra vertices
 // and a code indicating how the split should be made.
 p0=(long)(v0->ip[index]);
 p1=(long)(v1->ip[index]);
 p2=(long)(v2->ip[index]);
 if(IN_LESS_EQUAL_VOXEL(p0,p) && IN_LESS_EQUAL_VOXEL(p1,p) && IN_LESS_EQUAL_VOXEL(p2,p)){
   *vxi = 1; // code for less voxel
   return 0; // no extra vertices required
 }
 else if(IN_GREATER_EQUAL_VOXEL(p0,p) && IN_GREATER_EQUAL_VOXEL(p1,p) && IN_GREATER_EQUAL_VOXEL(p2,p)){
   *vxi = 2; // code for greater voxel
   return 0; // no extra vertices required
 }
 else{ // face must cross the voxel boundary
   // these require 1 additional vertex and face
   if     (IN_PLANE(p0,p)){if(IN_LESS_VOXEL(p2,p))*vxi = 3; else *vxi = 13;  return 1;}
   else if(IN_PLANE(p1,p)){if(IN_LESS_VOXEL(p0,p))*vxi = 4; else *vxi = 14;  return 1;}
   else if(IN_PLANE(p2,p)){if(IN_LESS_VOXEL(p1,p))*vxi = 5; else *vxi = 15;  return 1;}
   // these require 2 additional vertices and faces
   else if(IN_LESS_VOXEL(p0,p)){
     if     (IN_LESS_VOXEL(p1,p)){*vxi = 12; return 2;} 
     else if(IN_LESS_VOXEL(p2,p)){*vxi = 11; return 2;} 
     else                        {*vxi = 6;  return 2;}  
   }
   else if(IN_LESS_VOXEL(p1,p)){
     if     (IN_LESS_VOXEL(p2,p)){*vxi = 10; return 2;} 
     else if(IN_LESS_VOXEL(p0,p)){*vxi = 12; return 2;} 
     else                        {*vxi = 7;  return 2;}  
   }
   else if(IN_LESS_VOXEL(p2,p)){
     if     (IN_LESS_VOXEL(p0,p)){*vxi = 11; return 2;} 
     else if(IN_LESS_VOXEL(p1,p)){*vxi = 10; return 2;} 
     else                        {*vxi = 8;  return 2;}  
   }
 }
 //if(debug != NULL)fprintf(debug,"Voxel Split - this point should not be reached\n");
 * vxi = 0; 
 return 0;
}

static BOOL SplitVoxel(node *np, node *na, node *nb, long p, long index){ // index = 2 for z coord 0 for x coord
  // split polygons in node "np" into the two sub nodes "na" and "nb" at plane "p" 
 long i,nnp,nx,n1,n2,o,f,c,sz,xvf,nl,nr,ns,ng;
 tface *F;
 tvertex *V,*v0,*v1,*v2;
 nnp=np->n;
 // assign max and min for the new voxels
 //if(debug != NULL)fprintf(debug,"Split voxel with %ld faces in plane %ld at %ld\n",nnp,index,p); 
 for(i=0;i<3;i++){
   na->mi[i]=nb->mi[i]=np->mi[i];
   na->ma[i]=nb->ma[i]=np->ma[i];
 }
 na->ma[index]=p;     //  "na" points at voxel with the lesser coord value  
 nb->mi[index]=p;     //  "nb" to one with greater coord value
 // calculate number of extra vertices and faces required  
 for(i=0;i<ObjectCount;i++)ObjectXlist[i]=0;
 for(i=0;i<nnp;i++){  
   o = (long)(*(np->oid+i));
   f = (long)(*(np->fid+i));
   F=(Object[o].Ftracebase+f);
   V=Object[o].Vtracebase;
   v0=(V + F->V[0]);
   v1=(V + F->V[1]);
   v2=(V + F->V[2]);
   // calculate the number of extra vertices created due to possible face split and where they should go
   ObjectXlist[o] += FaceIntersectVoxelPlane(v0,v1,v2,p,index,(np->vxi+i)); // increment extra face/vertex count for object "o"
   //if(debug != NULL)fprintf(debug,"Object %ld Face %ld  Icode=%ld - vertices %ld %ld %ld \n",o,i,(long)(*(np->vxi+i)),F->V[0],F->V[1],F->V[2]);
 }
 //if(debug != NULL)for(i=0;i<ObjectCount;i++)fprintf(debug,"Object %ld has %ld extra faces in voxel\n",i,ObjectXlist[i]); 
 // allocate space for extra faces and vertices - this is absolute
 for(i=0,nx=0;i<ObjectCount;i++)if((xvf=ObjectXlist[i]) > 0){
   sz=(Object[i].NoTraceVertices+xvf)*(long)sizeof(tvertex);
   if((V=(tvertex  *)X__Realloc(Object[i].Vtracebase,sz)) == NULL)goto FAILED;
   Object[i].Vtracebase=V;
   sz=(Object[i].NoTraceFaces+xvf)*(long)sizeof(tface);
   if((F=(tface  *)X__Realloc(Object[i].Ftracebase,sz)) == NULL)goto FAILED;
   Object[i].Ftracebase=F;
   nx += xvf;
 }
 // allocate space for list of vertices and in child nodes. This is relative to voxel.
 // this will be over estimate since some of the added faces will go in one voxel and
 // some in the other
 ng=nnp+nx;
 //fprintf(debug,"Guessed voxel count %ld\n",ng);
 if(ng > 0){
   if((na->oid=(UC *)X__Malloc(sizeof(UC)*ng)) == NULL)goto FAILED;
   if((na->fid=(US *)X__Malloc(sizeof(US)*ng)) == NULL)goto FAILED;
   if((na->vxi=(UC *)X__Malloc(sizeof(UC)*ng)) == NULL)goto FAILED;
   if((nb->oid=(UC *)X__Malloc(sizeof(UC)*ng)) == NULL)goto FAILED;
   if((nb->fid=(US *)X__Malloc(sizeof(US)*ng)) == NULL)goto FAILED;
   if((nb->vxi=(UC *)X__Malloc(sizeof(UC)*ng)) == NULL)goto FAILED;
for(i=0;i<ng;i++){na->oid[i]=999,na->fid[i]=999,na->vxi[i]=999;nb->oid[i]=666,nb->fid[i]=666,nb->vxi[i]=666;} //remove when not debugging
 }
 // split vertices and faces - allocate them to the appropriate sub voxel - use VXI code
 nl=nr=ns=0;
 for(i=0,n1=0,n2=0;i<nnp;i++){  
   c = (long)(*(np->vxi+i));   
   o = (long)(*(np->oid+i));
   f = (long)(*(np->fid+i));
   if(c == 1){  // assigned to lower voxel
     *(na->fid+n1)=f; 
     *(na->oid+n1)=o; n1++;
     nl++;
   }
   else if(c == 2){
     *(nb->fid+n2)=f; 
     *(nb->oid+n2)=o; n2++;
     nr++;
   }
   else{ // split face
     F=(Object[o].Ftracebase+f);
     V=Object[o].Vtracebase;
     v0=(V + F->V[0]);
     v1=(V + F->V[1]);
     v2=(V + F->V[2]);
     SplitFaceAtVoxelPlane(o,f,c,V,Object[o].Ftracebase,v0,v1,v2,F,p,index,na->oid,na->fid,nb->oid,nb->fid,&n1,&n2);
     ns++;
   }
 }
 //fprintf(debug,"Assigned to lower %ld  to upper %ld  and split %ld\n",nl,nr,ns);
 // assign the correct sizes to voxel face list
 if(ng > 0 && n1 == 0){
   X__Free(na->oid); na->oid=NULL;
   X__Free(na->fid); na->fid=NULL;
   X__Free(na->vxi); na->vxi=NULL;
 }
 else if(n1 > 0){
   na->oid=(UC *)X__Realloc(na->oid,sizeof(UC)*n1);
   na->vxi=(UC *)X__Realloc(na->vxi,sizeof(UC)*n1);
   na->fid=(US *)X__Realloc(na->fid,sizeof(US)*n1);
 }
 if(ng > 0 && n2 == 0){
   X__Free(nb->oid); nb->oid=NULL;
   X__Free(nb->fid); nb->fid=NULL;
   X__Free(nb->vxi); nb->vxi=NULL;
 }
 if(n2 > 0){
   nb->oid=(UC *)X__Realloc(nb->oid,sizeof(UC)*n2);
   nb->vxi=(UC *)X__Realloc(nb->vxi,sizeof(UC)*n2);
   nb->fid=(US *)X__Realloc(nb->fid,sizeof(US)*n2);
 }
 // parent voxel now split so free face list
 if(np->oid != NULL)X__Free(np->oid); np->oid=NULL;
 if(np->fid != NULL)X__Free(np->fid); np->fid=NULL;
 if(np->vxi != NULL)X__Free(np->vxi); np->vxi=NULL;
 na->n=n1;
 nb->n=n2;
 //fprintf(debug,"%ld split into %ld  %ld\n",np->n, na->n,nb->n); fflush(debug); 
 np->n=0; 
 return OK;
 FAILED:
 return FAIL;
}

static int MakeVoxels(node *np, int level){  // main recursive function
 static long i,j,k,o,f,n;
 static double x,y,z,an;
 static tvertex *v,*v0,*v1,*v2;
 static tface   *F;
 static vector p;
 static node Nd[6]; // used for temporary datastructures
 RenderYield();
 if(R_terminator() == FAIL)longjmp(x_buf,-1);
 //if(debug != NULL)fprintf(debug,"Examining voxel %ld with %ld faces - depth=%ld\n",np->id,np->n,tree_depth); 
 if(level > threshold_r)return OK;
 if(np->n < threshold_f)return OK;   /* reached voxel minimun mumber of faces */
 if((np->ma[0]-np->mi[0]) < threshold_z ||
    (np->ma[1]-np->mi[1]) < threshold_z ||
    (np->ma[2]-np->mi[2]) < threshold_z)return OK;  // minimum size of voxels
 //if(debug != NULL){
 //fprintf(debug,"Splitting!  bounds [%ld %ld] [%ld %ld] [%ld %ld] size = %ld %ld %ld\n",
 //np->mi[0],np->ma[0],np->mi[1],np->ma[1],np->mi[2],np->ma[2],
 //np->ma[0]-np->mi[0],np->ma[1]-np->mi[1],np->ma[2]-np->mi[2]);
 //if(np->last != NULL)fprintf(debug,"Parent %ld\n",np->last->id); else fprintf(debug,"Spliting ROOT\n");
 //}
#if 0
 x=(double)(np->mi[0]+np->ma[0])*0.5;
 y=(double)(np->mi[1]+np->ma[1])*0.5;
 z=(double)(np->mi[2]+np->ma[2])*0.5;
#else
 {// This attempts to place the dividing plane through a vertex as close to the center as
  // possible. THis avoids subdividing some faces.
  // It assumes that every face lies in the voxel and all vertices are either
  // on the edge or inside. 
  long ds,ip,md[3],mxx[3]={BIGP,BIGP,BIGP};
  tvertex *vid[3]={NULL,NULL,NULL};
  md[0]=(np->mi[0]+np->ma[0])/2;
  md[1]=(np->mi[1]+np->ma[1])/2;
  md[2]=(np->mi[2]+np->ma[2])/2;
  for(i=0;i<np->n;i++){ // each face
    o = (long)(*(np->oid+i));   f = (long)(*(np->fid+i));
    F=(Object[o].Ftracebase+f); v=Object[o].Vtracebase;
    for(j=0;j<3;j++){ // each vertex
      v0=(v + F->V[j]);
      for(k=0;k<3;k++){//each coord 
        ip=v0->ip[k]; 
        if(ip > np->mi[k]+16 && ip < np->ma[k]-16 & (ds=abs(ip - md[k])) < mxx[k]){ // away from box wall
          mxx[k]=ds; vid[k]=v0; // a suitable vertex
        }
      }
    }  
  } // all faces
  for(k=0;k<3;k++){ // check to see if this vertex is suitable for split or just use the mid point 
    if(vid[k] != NULL){
     ip=vid[k]->ip[k];
     if(ip > np->mi[k]+128 && ip < np->ma[k]-128)mxx[k]=ip; else mxx[k]=md[k];
    } 
  }
  x=(double)mxx[0]; y=(double)mxx[1]; z=(double)mxx[2]; 
 }
#endif
 // x,y,z are the coordinates of the voxel splitting plane 

 // Partition the faces in Voxel "np"  by spltting at the split plane
 // first allocating the 8 child nodes in this octree.
 if(level > level_depth_max)level_depth_max=level;
 // temporary octree structures
 for(i=0;i<6;i++){
   Nd[i].n=0;
   Nd[i].oid=NULL;
   Nd[i].fid=NULL;
   Nd[i].vxi=NULL;
 }
 //  Create the child nodes
 level += 1;
 for(i=0;i<8;i++){
   if((np->octree[i]=(node *)X__Malloc(sizeof(node))) == NULL)return FAIL;
   np->octree[i]->n=0;
   np->octree[i]->oid=NULL;
   np->octree[i]->fid=NULL;
   np->octree[i]->vxi=NULL;
   for(j=0;j<8;j++){
     np->octree[i]->octree[j]=NULL;
   }
   kcount++;
   np->octree[i]->id=kcount;   // used for testing // RSFd
   np->octree[i]->last=np;     // testing
   np->octree[i]->level=level; // testing
 }
 if(trace_partition == 0){
 // split along planes - adding new faces this will re-create the vertex & face lists
 if(!SplitVoxel(np, &Nd[0], &Nd[1],(long)z,2))return FAIL; // give up if split failed
 if(!SplitVoxel(&Nd[0],&Nd[2],&Nd[3],(long)x,0))return FAIL;
 if(!SplitVoxel(&Nd[1],&Nd[4],&Nd[5],(long)x,0))return FAIL;
 if(!SplitVoxel(&Nd[2],np->octree[0],np->octree[1],(long)y,1))return FAIL;
 if(!SplitVoxel(&Nd[3],np->octree[2],np->octree[3],(long)y,1))return FAIL;
 if(!SplitVoxel(&Nd[4],np->octree[4],np->octree[5],(long)y,1))return FAIL;
 if(!SplitVoxel(&Nd[5],np->octree[6],np->octree[7],(long)y,1))return FAIL;
 } else {
 if(!PartitionVoxel(np, &Nd[0], &Nd[1],(long)z,2))return FAIL; // give up if split failed
 if(!PartitionVoxel(&Nd[0],&Nd[2],&Nd[3],(long)x,0))return FAIL;
 if(!PartitionVoxel(&Nd[1],&Nd[4],&Nd[5],(long)x,0))return FAIL;
 if(!PartitionVoxel(&Nd[2],np->octree[0],np->octree[1],(long)y,1))return FAIL;
 if(!PartitionVoxel(&Nd[3],np->octree[2],np->octree[3],(long)y,1))return FAIL;
 if(!PartitionVoxel(&Nd[4],np->octree[4],np->octree[5],(long)y,1))return FAIL;
 if(!PartitionVoxel(&Nd[5],np->octree[6],np->octree[7],(long)y,1))return FAIL;
 }
 /* recurse on to next level */
 MakeVoxels(np->octree[0],level);
 MakeVoxels(np->octree[1],level);
 MakeVoxels(np->octree[2],level);
 MakeVoxels(np->octree[3],level);
 MakeVoxels(np->octree[4],level);
 MakeVoxels(np->octree[5],level);
 MakeVoxels(np->octree[6],level);
 MakeVoxels(np->octree[7],level);
 return OK;
}

int buildBST(void){
 char temp[256];
 long i,j,k,l,coord,nv=0;
 long divisions;
 tvertex *v;
 tface   *f;
 nfaces=0; kcount=0; level_depth_max=0;
 voxel_xmin=voxel_ymin=voxel_zmin=BIGP;
 voxel_xmax=voxel_ymax=voxel_zmax=BIGN;
 Render_Message(31,0,NULL);
 ROOT=NULL;
 if(ObjectCount > 0)for(i=0;i<ObjectCount;i++){
   if(!Object[i].in_use)continue;
   Object[i].NoTraceVertices=Object[i].NoVertices;
   Object[i].NoTraceFaces=Object[i].NoFaces;
   nfaces += Object[i].NoTraceFaces;
   nv += Object[i].NoTraceVertices;
   // copy vertex and face data to the geometry data structures for ray-tracing
   if((Object[i].Vtracebase=(tvertex *)X__Malloc(Object[i].NoTraceVertices*sizeof(tvertex))) == NULL)return FAIL;
   for(j=0;j<Object[i].NoTraceVertices;j++){
     VECCOPY(Object[i].Vbase[j].ip,Object[i].Vtracebase[j].ip)
   }
   if((Object[i].Ftracebase=(tface *)X__Malloc(Object[i].NoTraceFaces*sizeof(tface))) == NULL)return FAIL;
   for(j=0;j<Object[i].NoTraceFaces;j++){
     VECCOPY(Object[i].Fbase[j].V,Object[i].Ftracebase[j].V)
     Object[i].Ftracebase[j].o=i;
     Object[i].Ftracebase[j].f=j;
     { // calculate normal vector used in partition algorithm
     vector v1,v2,n;
     VECSUB(Object[i].Vtracebase[Object[i].Ftracebase[j].V[1]].ip,
            Object[i].Vtracebase[Object[i].Ftracebase[j].V[0]].ip,v1)
     VECSUB(Object[i].Vtracebase[Object[i].Ftracebase[j].V[2]].ip,
            Object[i].Vtracebase[Object[i].Ftracebase[j].V[0]].ip,v2)
     CROSS(v1,v2,n)
     normalize(n);
     VECCOPY(n,Object[i].Ftracebase[j].n)
     }
   }
 }
 if(nfaces == 0)return 0;
 if((ROOT=(node *)X__Malloc(sizeof(node))) == NULL)return FAIL;
 ROOT->n=0;
 if((ROOT->oid=(UC *)X__Malloc(sizeof(UC ) * nfaces)) == NULL){
   destroyBST(ROOT); return FAIL;
 }
 if((ROOT->fid=(US *)X__Malloc(sizeof(US ) * nfaces)) == NULL){
   destroyBST(ROOT); return FAIL;
 }
 if((ROOT->vxi=(UC *)X__Malloc(sizeof(UC ) * nfaces)) == NULL){
   destroyBST(ROOT); return FAIL;
 }
 ROOT->n=nfaces;  // all polygons in root node to start with
 if((ObjectXlist=(long *)X__Malloc(sizeof(long)*ObjectCount)) == NULL)return FAIL;
 for(i=0,k=0;i<ObjectCount;i++){
  if(!Object[i].in_use)continue;
  f=Object[i].Ftracebase;
   for(j=0;j<Object[i].NoTraceFaces;j++){
     *(ROOT->oid+k)=(UC)i;
     *(ROOT->fid+k)=(US)j;
     *(ROOT->vxi+k)=(UC)0;
     k++;
     f++;
   }
 }
 for(i=0;i<8;i++)ROOT->octree[i]=NULL;
 ROOT->mi[0]=ROOT->mi[1]=ROOT->mi[2]=BIGP;
 ROOT->ma[0]=ROOT->ma[1]=ROOT->ma[2]=BIGN;
 ROOT->id=kcount; ROOT->last=NULL;  // Used for testing //RSFd
 for(i=0,k=0;i<ObjectCount;i++){
   if(!Object[i].in_use)continue;
   v = Object[i].Vtracebase;
   for(j=0;j<Object[i].NoTraceVertices;j++){
     for(l=0;l<3;l++){  //p are the transformed coordinates after effects and clipping
       if((coord=(long)(v->ip[l])) < ROOT->mi[l])ROOT->mi[l]=coord;
       if( coord                   > ROOT->ma[l])ROOT->ma[l]=coord;
     }
     k++;
     v++;
   }
 }
 if(k == 0)return 0;
 for(i=0;i<3;i++){
   ROOT->mi[i] -= 2;  // to ensure voxel contains all vertices
   ROOT->ma[i] += 2;
 }
 longest_ray = max((ROOT->ma[1]-ROOT->mi[1]),(ROOT->ma[0]-ROOT->mi[0]));
 longest_ray = max(longest_ray,(ROOT->ma[2]-ROOT->mi[2]));
 threshold_r=trace_depth;             //THRESHOLDR;
 threshold_f=trace_face_count_target; //THRESHOLDF;
 divisions = (long)(pow(2.0,(double)threshold_r));  // max number of subdivisions along a side;
 threshold_z=(long)(longest_ray/divisions);
 longest_ray *= 1.4142; 
 //if(debug != NULL){ 
 //  fprintf(debug,"Dimensions x %ld  y %ld  z%ld\nLongest Ray %lf\n",ROOT->ma[0]-ROOT->mi[0],
 //        ROOT->ma[1]-ROOT->mi[1],ROOT->ma[2]-ROOT->mi[2],longest_ray);
 //  fprintf(debug,"Thr_r %d  Thr_f %d Thr_d  %d  (div %d)\n",threshold_r,threshold_f,threshold_z,divisions);
 //}
 /* make the voxels   --  at this time ALL the faces are in voxel ROOT */
 if(setjmp(x_buf) != 0 || MakeVoxels(ROOT,0) == FAIL){
   destroyBST(ROOT);
   X__Free(ObjectXlist);
   return FAIL;
 }
 X__Free(ObjectXlist);
 i=nfaces;
 nfaces=0;
 maxfaces=0;
 minfaces=65536*65536;
 eliminateDuplicatesInBST(ROOT);
 countBST(ROOT);
 //printBST(ROOT);  
 // SaveDivided(nvoxels,ROOT,0,"c:\\divided_model.mfx"); // RSFd
 sprintf(temp,
  "Scene: F=%ld V=%ld [Dp=%ld Th=%ld F/v=%ld] Size %ld(%ld) [Vox=%ld Fc=%ld]",
        i,nv,level_depth_max,threshold_f,maxfaces,
        max(min(voxel_xmax,voxel_ymax),voxel_zmax)/512,
        min(min(voxel_xmin,voxel_ymin),voxel_zmin)/512,
        kcount,nfaces);
 WindowedRenderInfo(temp);
 //if(debug != NULL){ 
 //  fprintf(debug,"Number of original faces %d  Number of original vertices %d\n",i,nv);
 //  fprintf(debug,"Number of faces in structured model %d\n",nfaces);
 //  fprintf(debug,"Number of Voxels %d\n",kcount);
 //  fprintf(debug,"Max number of faces in Voxel %d\n",maxfaces); 
 //}
 return 0;
}

void freeBST(void){
 int i;
 if(ROOT != NULL)destroyBST(ROOT);
 if(ObjectCount > 0)for(i=0;i<ObjectCount;i++){
   if(!Object[i].in_use)continue;
   if(Object[i].Vtracebase != NULL)X__Free(Object[i].Vtracebase); 
   Object[i].Vtracebase=NULL;  Object[i].NoTraceVertices=0;
   if(Object[i].Ftracebase != NULL)X__Free(Object[i].Ftracebase); 
   Object[i].Ftracebase=NULL; Object[i].NoTraceFaces=0;
 }
}

static void destroyBST(node *np){
 if(np->octree[0] != NULL){destroyBST(np->octree[0]); np->octree[0]=NULL;}
 if(np->octree[1] != NULL){destroyBST(np->octree[1]); np->octree[1]=NULL;}
 if(np->octree[2] != NULL){destroyBST(np->octree[2]); np->octree[2]=NULL;}
 if(np->octree[3] != NULL){destroyBST(np->octree[3]); np->octree[3]=NULL;}
 if(np->octree[4] != NULL){destroyBST(np->octree[4]); np->octree[4]=NULL;}
 if(np->octree[5] != NULL){destroyBST(np->octree[5]); np->octree[5]=NULL;}
 if(np->octree[6] != NULL){destroyBST(np->octree[6]); np->octree[6]=NULL;}
 if(np->octree[7] != NULL){destroyBST(np->octree[7]); np->octree[7]=NULL;}
 if(np->fid != NULL){X__Free(np->fid); np->fid=NULL;}
 if(np->oid != NULL){X__Free(np->oid); np->oid=NULL;}
 if(np->vxi != NULL){X__Free(np->vxi); np->vxi=NULL;}
 X__Free(np); np=NULL;
 return;
}

static void countBST(node *np){
 char temp[256];
 int i;
 if(np->ma[0]-np->mi[0] < voxel_xmin)voxel_xmin=np->ma[0]-np->mi[0];
 if(np->ma[1]-np->mi[1] < voxel_ymin)voxel_ymin=np->ma[1]-np->mi[1];
 if(np->ma[2]-np->mi[2] < voxel_zmin)voxel_zmin=np->ma[2]-np->mi[2];
 if(np->ma[0]-np->mi[0] > voxel_xmax)voxel_xmax=np->ma[0]-np->mi[0];
 if(np->ma[1]-np->mi[1] > voxel_ymax)voxel_ymax=np->ma[1]-np->mi[1];
 if(np->ma[2]-np->mi[2] > voxel_zmax)voxel_zmax=np->ma[2]-np->mi[2];
 if(np->octree[0] != NULL)countBST(np->octree[0]);
 if(np->octree[1] != NULL)countBST(np->octree[1]);
 if(np->octree[2] != NULL)countBST(np->octree[2]);
 if(np->octree[3] != NULL)countBST(np->octree[3]);
 if(np->octree[4] != NULL)countBST(np->octree[4]);
 if(np->octree[5] != NULL)countBST(np->octree[5]);
 if(np->octree[6] != NULL)countBST(np->octree[6]);
 if(np->octree[7] != NULL)countBST(np->octree[7]);
 if(np->fid != NULL){
   if(np->n > maxfaces)maxfaces=np->n;
   if(np->n < minfaces)minfaces=np->n;
   nfaces += np->n;
 }
 return;
}

static void eliminateDuplicatesInBST(node *n){ // this does not reduce the number of divided Trace faces 
 if(n->n > 0){                                 // it just ensures that the list of each trace face in a voxel does not 
   static long i,o,f,ft,pfid,poid,c;           // point at same model face more than once. Helps speed up the tracing.
   US *next;
   pfid=-1,poid=-1; c=0;
   if((next=(US *)X__Malloc(sizeof(US)*n->n)) == NULL)return;
   for(i=0;i<n->n;i++){
     o=n->oid[i];
     f=n->fid[i];
     ft=(Object[o].Ftracebase + f)->f;  // points at equivalent face in main data structures
     if(!(ft == pfid && o == poid))next[c++]=i;  // lists first of any runs !
     pfid=ft; poid=o;
   }
   if(c < n->n)for(i=0;i<c;i++){ // re-order the list remove duplicates and shorten
     n->oid[i]=n->oid[next[i]];  // next[i] is always >= i
     n->fid[i]=n->fid[next[i]];
   }
   n->n=c;
   X__Free(next); 
 }
 if(n->octree[0] != NULL)eliminateDuplicatesInBST(n->octree[0]);
 if(n->octree[1] != NULL)eliminateDuplicatesInBST(n->octree[1]);
 if(n->octree[2] != NULL)eliminateDuplicatesInBST(n->octree[2]);
 if(n->octree[3] != NULL)eliminateDuplicatesInBST(n->octree[3]);
 if(n->octree[4] != NULL)eliminateDuplicatesInBST(n->octree[4]);
 if(n->octree[5] != NULL)eliminateDuplicatesInBST(n->octree[5]);
 if(n->octree[6] != NULL)eliminateDuplicatesInBST(n->octree[6]);
 if(n->octree[7] != NULL)eliminateDuplicatesInBST(n->octree[7]);
 return;

}

#define TOL    0.00001     

// This group of functions tracks a ray by working down the voxel tree recursively and finding which voxels
// it intersects. When we find that the ray passes through a voxel with some faces we check them for intersection
// with the whole ray. -  Seems to work with shadows 

static BOOL RayIntersectVoxelSides(node *n, double ddn, vector p1, vector d, vector p2, long i, long j, long k){
 double ppn,u,v,mu;
 if((p1[i] <= n->mi[i] && p2[i] >= n->mi[i]) ||
    (p2[i] <= n->mi[i] && p1[i] >= n->mi[i])){
   ppn=n->mi[i]-p1[i]; mu=ppn/ddn;
   u=p1[j]+mu*d[j];
   v=p1[k]+mu*d[k];
   if(u >= n->mi[j] && u <= n->ma[j] && v >= n->mi[k] && v <= n->ma[k])return TRUE;
 }
 if((p1[i] <= n->ma[i] && p2[i] >= n->ma[i]) ||
    (p2[i] <= n->ma[i] && p1[i] >= n->ma[i])){
   ppn=n->ma[i]-p1[i]; mu=ppn/ddn;
   u=p1[j]+mu*d[j];
   v=p1[k]+mu*d[k];
   if(u >= n->mi[j] && u <= n->ma[j] && v >= n->mi[k] && v <= n->ma[k])return TRUE;
 }
 return FALSE;
}

static BOOL RayThroughVoxel(node *n, vector p1, vector d, vector p2){
 double ddn,ppn,u,v;
 if(p1[0] < n->mi[0] && p2[0] < n->mi[0])return FALSE;
 if(p1[1] < n->mi[1] && p2[1] < n->mi[1])return FALSE;
 if(p1[2] < n->mi[2] && p2[2] < n->mi[2])return FALSE;
 if(p1[0] > n->ma[0] && p2[0] > n->ma[0])return FALSE;
 if(p1[1] > n->ma[1] && p2[1] > n->ma[1])return FALSE;
 if(p1[2] > n->ma[2] && p2[2] > n->ma[2])return FALSE;
 if(p1[0] >= n->mi[0] && p1[0] <= n->ma[0] &&
    p1[1] >= n->mi[1] && p1[1] <= n->ma[1] &&
    p1[2] >= n->mi[2] && p1[2] <= n->ma[2])return TRUE;
 if(p2[0] >= n->mi[0] && p2[0] <= n->ma[0] &&
    p2[1] >= n->mi[1] && p2[1] <= n->ma[1] &&
    p2[2] >= n->mi[2] && p2[2] <= n->ma[2])return TRUE;
 ddn=d[0]; if(fabs(ddn) > TOL)if(RayIntersectVoxelSides(n,ddn,p1,d,p2,0,1,2))return TRUE;
 ddn=d[1]; if(fabs(ddn) > TOL)if(RayIntersectVoxelSides(n,ddn,p1,d,p2,1,2,0))return TRUE;
 ddn=d[2]; if(fabs(ddn) > TOL)if(RayIntersectVoxelSides(n,ddn,p1,d,p2,2,0,1))return TRUE;
 return FALSE;
}

static void FeelerHitsFace(node *np, vector p, vector d, vector pl, double ld, 
                           long obj, long *glassc, BOOL *hit, BOOL *abort){
 int i,j,o;
 face   *f;
 long  ft;
 vertex *v;
 matl   *mat;
 double depth;
 unsigned char  mattr;
 long glass = 0;
 node *nx;
 if(*abort)return;
 //if(np->last != NULL) fprintf(debug,"Testing voxel %ld parent %ld (%ld faces)\n",np->id,np->last->id,np->n);
 //else fprintf(debug,"Testing root voxel %ld \n",np->id);
 if(!RayThroughVoxel(np,p,d,pl))return;
 //fprintf(debug,"Ray passes through voxel With %ld faces\n",np->n);
 if(np->n > 0){ /* some faces */
   for(i=0;i<np->n;i++){
     o=(int)(*(np->oid + i)); // change to reflect the actual face (this should be same for both main and trace lists
     if(!Object[o].cast_shadow)continue;
     if(!Object[obj].self_shadow && obj == o)continue;
     ft=(Object[o].Ftracebase + (int)(*(np->fid + i)))->f;  // points at equivalent face in the original model
     f=Object[o].Fbase + (int)(ft);   // change to relect the actual face
     v=Object[o].Vbase;
     if(f->ffmat >= 0 && f->ffmat < Object[o].NoMats){
       if(Object[o].Tbase != NULL){
         mat=(Object[o].Tbase + f->ffmat);
         mattr=mat->transp;
       }
       else mattr=1;
     }
     else mattr=1;
     //fprintf(debug,"Testing face %ld object %ld face %ld\n",i,o,(int)(*(np->fid + i)));
     if(R_intersect_face(p,d,
              (v + f->V[0])->p,(v + f->V[1])->p,(v + f->V[2])->p,
               f->n,&depth,32.0) && (depth < ld)){ // need this so not on far side of light
       if(mattr > 1)glassc++;
       else{
         *abort=TRUE;
         *hit=TRUE;
         return;
       }           
     }
   }
 }
 if((nx=np->octree[0]) != NULL)FeelerHitsFace(nx,p,d,pl,ld,obj,glassc,hit,abort);  else return; // all the rest will be NULL if this one is
 if((nx=np->octree[1]) != NULL)FeelerHitsFace(nx,p,d,pl,ld,obj,glassc,hit,abort); 
 if((nx=np->octree[2]) != NULL)FeelerHitsFace(nx,p,d,pl,ld,obj,glassc,hit,abort); 
 if((nx=np->octree[3]) != NULL)FeelerHitsFace(nx,p,d,pl,ld,obj,glassc,hit,abort);
 if((nx=np->octree[4]) != NULL)FeelerHitsFace(nx,p,d,pl,ld,obj,glassc,hit,abort);
 if((nx=np->octree[5]) != NULL)FeelerHitsFace(nx,p,d,pl,ld,obj,glassc,hit,abort);
 if((nx=np->octree[6]) != NULL)FeelerHitsFace(nx,p,d,pl,ld,obj,glassc,hit,abort);
 if((nx=np->octree[7]) != NULL)FeelerHitsFace(nx,p,d,pl,ld,obj,glassc,hit,abort);
 return;
}

short check_for_shadow(long obj, vector p, vector natp){
/* p is point for shadow   natp  is normal at point p */
 int i,j;
 vector d,pp;
 double ld,alpha;
 long glassc=0;
 BOOL abort,hit;

 for(j=0;j<Nlights;j++)if(Lights[j].type != NOSHADOW &&
                          Lights[j].type != DUMMYL){
   if(Lights[j].type == SPOTLIGHT && Lights[j].shadow_buffer != NULL)continue;
   VECSUB(Lights[j].p,p,d)
   if(DOT(p,natp)*DOT(d,natp) > 0)continue; /* not illuminated by light */
   ld=(double)((d[0]*d[0]) + (d[1]*d[1]) + (d[2]*d[2]));
   if(ld > 1.e-10){
     ld=sqrt(ld);
     VECSCALE((1.0/ld),d,d)
     if(Lights[j].type == SPOTLIGHT){
       if((alpha=DOT(d,Lights[j].d)) < Lights[j].dot2){
            continue; /* outside light cone */
       }
     }
     abort=FALSE;
     hit=FALSE;
     VECCOPY(p,pp)     
     //fprintf(debug,"Starting test from point (%lf %lf %lf) to light at (%lf %lf %lf)\n",p[0],p[1],p[2],Lights[j].p[0],Lights[j].p[1],Lights[j].p[2]);
     FeelerHitsFace(ROOT,pp,d,Lights[j].p,ld,obj,&glassc,&hit,&abort);
     if(hit)return -1;
   }
 }
 return glassc;
}

static void GetRayIntersectionFace(node *np, vector p, vector d, // start at p and in direction d
                           double *ld, // start off with longest distance to be sure outside object
                           long *closest_o, long *closest_f){
 // work with whole ray - work through whole of Octree structure - find any faces the ray intersects
 // only look for any faces that lie closer to the starting point when testing the next voxel level
 int i,j,o,fi;
 face   *f;
 long ft;
 vertex *v;
 double depth;
 node *nx;
 vector pl;
 VECSUM(p,(*ld)*d,pl)  // point on outside - end of ray 
 if(!RayThroughVoxel(np,p,d,pl))return;  // between p and pl
 if(np->n > 0){ /* some faces */
   for(i=0;i<np->n;i++){
     o=(int)(*(np->oid + i)); // change to reflect the actual face -- this should be the same - because no new objects are introduced by the face splitting
     ft=(Object[o].Ftracebase + (int)(*(np->fid + i)))->f;  // points at equivalent face in original model
     f=Object[o].Fbase + (fi=(int)(ft));   // change to relect the actual face
     v=Object[o].Vbase;
     if(R_intersect_face(p,d,
              (v + f->V[0])->p,(v + f->V[1])->p,(v + f->V[2])->p,
               f->n,&depth,32.0)){
       if(depth < *ld){
         *ld=depth;   // only check any subsequent voxels if they are nearer the origin
         *ld=depth - 10.0;    // only check any subsequent voxels if they are nearer the origin
         *closest_o = o;   // change to reflect  actual face ID
         *closest_f = fi;
       }           
     }
   }
 }
 if((nx=np->octree[0]) != NULL)GetRayIntersectionFace(nx,p,d,ld,closest_o,closest_f);
 if((nx=np->octree[1]) != NULL)GetRayIntersectionFace(nx,p,d,ld,closest_o,closest_f); 
 if((nx=np->octree[2]) != NULL)GetRayIntersectionFace(nx,p,d,ld,closest_o,closest_f); 
 if((nx=np->octree[3]) != NULL)GetRayIntersectionFace(nx,p,d,ld,closest_o,closest_f);
 if((nx=np->octree[4]) != NULL)GetRayIntersectionFace(nx,p,d,ld,closest_o,closest_f);
 if((nx=np->octree[5]) != NULL)GetRayIntersectionFace(nx,p,d,ld,closest_o,closest_f);
 if((nx=np->octree[6]) != NULL)GetRayIntersectionFace(nx,p,d,ld,closest_o,closest_f);
 if((nx=np->octree[7]) != NULL)GetRayIntersectionFace(nx,p,d,ld,closest_o,closest_f);
 return;
}

BOOL trace_starting_ray(vector d, vector p,
                        int *ob, int *fc, face **ff){
/*  p is point of ray origin on entry point of hit on exit
    d is incident direction of ray (normalized)
*/
 int o,q;
 double depth;
 o=-1; 
 depth=(double)FARAWAY;
 GetRayIntersectionFace(ROOT,p,d,&depth,&o,&q);
 if(o < 0)return FALSE; // no hit with this ray  
 *ff=(Object[o].Fbase + q);
 *fc=q;
 *ob=o;
 VECSCALE(depth,d,p)
 return TRUE;
}

short trace_reflection_ray(vector p, vector d, /* d is normalized */
                           double depth_max, double *colour,
                           int at_o, face *at_f,
                           long trace_reflection_depth,
                           long trace_refraction_depth,
                           long refractive_flag){
 int o=-1,q;
 face *f;
 vertex *v;
 double depth=longest_ray;
TRY_START_AGAIN:
 o=-1; 
 depth=longest_ray;
 VECSUM(p,128.0*d,p) // under test
 GetRayIntersectionFace(ROOT,p,d,&depth,&o,&q);
 if(o < 0)return 0; // no hit with this ray  
 f=(Object[o].Fbase + q);
 if(o == at_o && f == at_f){
   // hit the same face - increment the starting point along the ray and try again
   goto TRY_START_AGAIN; 
 }
 v=Object[o].Vbase;
 f=(Object[o].Fbase + q);
 depth_max=depth;
 if(trace_reflection_depth < trace_depth_max){
   //colour[0] = (double)0;  colour[1] = (double)255; colour[2] = (double)0; 
   trace_reflection_depth++;
   GetSurfaceValue(o,f,v,p,d,
   colour,(colour+1),(colour+2),
   trace_reflection_depth,
   trace_refraction_depth,
   refractive_flag);
 }
 return 1;
}

void trace_transmission_ray(vector p, vector d, vector n,
                          int at_o, face *at_f,
                          double depth_max, double *colour,
                          long trace_reflection_depth,
                          long trace_refraction_depth,
                          long refractive_flag){
/*  p is point on glass face of ray hit
    d is incident direction of ray (normalized)
    n is normal of face on which ray is incident
    depth_max is maximum depth to trace to
*/
 int o,q;
 vector pc;
 face   *f;
 vertex *v;
 double depth;
 TRY_START_AGAIN:
 o=-1;
 depth=longest_ray;
 // go through the current face by 1000 units to make sure we dont hit the same face
 VECSUM(p,128.0*d,p) 
 GetRayIntersectionFace(ROOT,p,d,&depth,&o,&q);
 if(o < 0)return; // no hit with this ray  
 f=(Object[o].Fbase + q);
 if(o == at_o && f == at_f){
   // hit the same face - increment the starting point along the ray and try again
   goto TRY_START_AGAIN; 
 }
 v=Object[o].Vbase;
 if(trace_refraction_depth < trace_depth_max){
  trace_refraction_depth++;
   GetSurfaceValue(o,f,v,p,d,
                   colour,(colour+1),(colour+2),
                   trace_reflection_depth,
                   trace_refraction_depth,
                   refractive_flag);
 }
 return;
}


void trace_refraction_ray(vector p, vector d, vector n,
                          int at_o, face *at_f,
                          double depth_max, double *colour, int refx,
                          long trace_reflection_depth,
                          long trace_refraction_depth,
                          long refractive_flag){
/*  p is point on glass face of ray hit
    d is incident direction of ray (normalized)
    n is normal of face on which ray is incident
    depth_max is maximum depth to trace to
*/
 int o,q;
 vector pc,temp,qq;
 face   *f;
 vertex *v;
 double depth,ddn;
 vector dp,gp;
 double pdotn,ddotn,mu,qnorm;
 unsigned char red,green,blue;
 // colour[0] = (double)255;  colour[1] = (double)0; colour[2] = (double)0;  return; 
 ddn=DOT(d,n);
 if(ddn > 0.0){
   ddn = -ddn; n[0] = -n[0]; n[1] = -n[1]; n[2] = -n[2];
 }
 VECSCALE(ddn,n,temp)
 VECSUB(d,temp,temp)
 if(refractive_flag == 0){
   VECSCALE(1000.0/(double)refx,temp,qq)
 }
 else{
   VECSCALE((double)refx/1000.0,temp,qq)
 }
 refractive_flag ^= 1;  /* each intersection is assumed air/glass */
 qnorm=qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2];
 if(qnorm > 1.0){
   VECSCALE(sqrt(qnorm-1.0),n,temp)
   VECSUM(qq,temp,d)
   ddn=d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
   ddn=sqrt(ddn);
   VECSCALE((1.0/ddn),d,d)  /* d is normalised new ray direction */
 }
 else{
   VECSCALE(sqrt(1.0-qnorm),n,temp)
   VECSUB(qq,temp,d)
   ddn=d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
   ddn=sqrt(ddn);
   VECSCALE((1.0/ddn),d,d)  /* d is normalised new ray direction */
 }
 /* get the ground or sky value for this ray in case no other Xn's */
 if(R_Nground > 0 && fabs(ddotn=dot(Ground.n,d)) > 1.e-3){
   VECSUB(Ground.p,p,dp)
   pdotn=DOT(Ground.n,dp);
   mu=pdotn/ddotn;
   if(mu > 0.0){
     vecscale(mu,d,gp);
     vecsum(p,gp,gp);
     GetGroundLight(gp,&colour[0],&colour[1],&colour[2]);
     depth_max = mu;
   }
   else goto NEXTBLOCK;
 }
 else{
   NEXTBLOCK:
   colour[0] = (double)Sky.Zenith[0];
   colour[1] = (double)Sky.Zenith[1];
   colour[2] = (double)Sky.Zenith[2];
   // this is NOT correct for GRADED but is a small problem
   if(Sky.type == SKYMAPPED){
     ddotn=d[2];
     if(fabs(ddotn) > 1.e-2){
       vecsub(Sky.map.p,p,dp);
       pdotn=dp[2];
       mu=pdotn/ddotn;
       vecscale(fabs(mu),d,gp);
       vecsum(p,gp,gp);
       MapFromProjection(&Sky.map,TILE,EOPAQUE,gp,&colour[0],&colour[1],&colour[2]);
     }
   }
   else if(Sky.type == SKYCUBE){
     MapEnvironmentInDirection(d,&colour[0],&colour[1],&colour[2],&Sky.map);
   }
 }
 TRY_START_AGAIN:
 o=-1; 
 depth=longest_ray;
 // go through the current face by 1000 units to make sure we dont hit the same face
 VECSUM(p,128.0*d,p) 
 GetRayIntersectionFace(ROOT,p,d,&depth,&o,&q);
 if(o < 0)return; // no hit with this ray  
 f=(Object[o].Fbase + q);
 if(o == at_o && f == at_f){
   // hit the same face - increment the starting point along the ray and try again
   goto TRY_START_AGAIN; 
 }
 v=Object[o].Vbase;
 if(trace_refraction_depth < trace_depth_max){
   trace_refraction_depth++;
   GetSurfaceValue(o,f,v,p,d,
                   colour,(colour+1),(colour+2),
                   trace_reflection_depth,
                   trace_refraction_depth,
                   refractive_flag);
 }
 return;
}


#if 0

static void printBST(node *n){
 if(debug == NULL)return;
 if(n->n > 0){  // only print voxels with faces
   static int i,o,f,ft;
   fprintf(debug,"Voxel id=%5ld (Has %5ld Faces) (At level %5ld) Edges[%8ld][%8ld][%8ld]",n->id,n->n, n->level, n->ma[0]-n->mi[0],n->ma[1]-n->mi[1],n->ma[2]-n->mi[2]);
   if(n->last != NULL)fprintf(debug," -> %ld\n",n->last->id);
   else               fprintf(debug," ROOT \n");
   for(i=0;i<n->n;i++){
     o=n->oid[i];
     f=n->fid[i];
     Object[o].Ftracebase[f].x = n->id;  // used for debugging idetifies which voxel this face belongs to
     ft=(Object[o].Ftracebase + f)->f;  // points at equivalent face in main data structures
     fprintf(debug, "V%5ld F%3ld  -  (F)%6ld (O=%d) Maps to face %6ld\n",n->id,i,f,o,ft);
   }
 }

 if(n->octree[0] != NULL)printBST(n->octree[0]);
 if(n->octree[1] != NULL)printBST(n->octree[1]);
 if(n->octree[2] != NULL)printBST(n->octree[2]);
 if(n->octree[3] != NULL)printBST(n->octree[3]);
 if(n->octree[4] != NULL)printBST(n->octree[4]);
 if(n->octree[5] != NULL)printBST(n->octree[5]);
 if(n->octree[6] != NULL)printBST(n->octree[6]);
 if(n->octree[7] != NULL)printBST(n->octree[7]);
 return;
}

#endif


static BOOL  EdgeCrossBoxFaceX(ivector vi1, ivector vi2, long p, 
                              long ymin, long ymax,  long zmin, long zmax){
 long d1,d2,vo1,vo2,imu;
 double mu;
 if((imu=vi2[0]-vi1[0]) == 0)return FALSE;
 d1=vi1[0]-p; d2=vi2[0]-p;
 if((double)d1*(double)d2 >= 0.0)return FALSE; // edge does not cross the face
 mu=((double)(p-vi1[0]))/(double)imu;
 vo1=vi1[1]+(long)(mu*((double)(vi2[1]-vi1[1])));
 vo2=vi1[2]+(long)(mu*((double)(vi2[2]-vi1[2])));
 if(ymin <= vo1 && vo1 <= ymax && zmin <= vo2 && vo2 <= zmax)return TRUE; 
 return FALSE;
}

static BOOL  EdgeCrossBoxFaceY(ivector vi1, ivector vi2, long p, 
                              long zmin, long zmax,  long xmin, long xmax){
 long d1,d2,vo1,vo2,imu;
 double mu;
 if((imu=vi2[1]-vi1[1]) == 0)return FALSE;
 d1=vi1[1]-p; d2=vi2[1]-p;
 if((double)d1*(double)d2 >= 0.0)return FALSE; // edge does not cross the face
 mu=((double)(p-vi1[1]))/(double)imu;
 vo1=vi1[2]+(long)(mu*((double)(vi2[2]-vi1[2])));
 vo2=vi1[0]+(long)(mu*((double)(vi2[0]-vi1[0])));
 if(zmin <= vo1 && vo1 <= zmax && xmin <= vo2 && vo2 <= xmax)return TRUE; 
 return FALSE;
}

static BOOL  EdgeCrossBoxFaceZ(ivector vi1, ivector vi2, long p, 
                              long xmin, long xmax,  long ymin, long ymax){
 long d1,d2,vo1,vo2,imu;
 double mu;
 if((imu=vi2[2]-vi1[2]) == 0)return FALSE;
 d1=vi1[2]-p; d2=vi2[2]-p;
 if((double)d1*(double)d2 >= 0.0)return FALSE; // edge does not cross the face
 mu=((double)(p-vi1[2]))/(double)imu;
 vo1=vi1[0]+(long)(mu*((double)(vi2[0]-vi1[0])));
 vo2=vi1[1]+(long)(mu*((double)(vi2[1]-vi1[1])));
 if(xmin <= vo1 && vo1 <= xmax && ymin <= vo2 && vo2 <= ymax)return TRUE; 
 return FALSE;
}

static BOOL OnFaceEdge(node *n, long id, ivector ip){
 if      (id == 0 &&
  ip[1]==n->mi[1] && n->mi[0]<=ip[0] && ip[0]<=n->ma[0] && n->mi[2]<=ip[2] && ip[2]<=n->ma[2])return TRUE; 
 else if (id == 1 &&
  ip[1]==n->ma[1] && n->mi[0]<=ip[0] && ip[0]<=n->ma[0] && n->mi[2]<=ip[2] && ip[2]<=n->ma[2])return TRUE; 
 else if (id == 2 &&
  ip[0]==n->mi[0] && n->mi[1]<=ip[1] && ip[1]<=n->ma[1] && n->mi[2]<=ip[2] && ip[2]<=n->ma[2])return TRUE; 
 else if (id == 3 &&
  ip[0]==n->ma[0] && n->mi[1]<=ip[1] && ip[1]<=n->ma[1] && n->mi[2]<=ip[2] && ip[2]<=n->ma[2])return TRUE; 
 else if (id == 4 &&
  ip[2]==n->mi[2] && n->mi[0]<=ip[0] && ip[0]<=n->ma[0] && n->mi[1]<=ip[1] && ip[1]<=n->ma[1])return TRUE; 
 else if (id == 5 &&
  ip[2]==n->ma[2] && n->mi[0]<=ip[0] && ip[0]<=n->ma[0] && n->mi[1]<=ip[1] && ip[1]<=n->ma[1])return TRUE; 
 return FALSE;
}

static BOOL BoxEdgeCrossFace(ivector p1, ivector p2, long index, long len, vector n, 
                             ivector pf0, ivector pf1, ivector pf2){
 vector d,v1,pi,w,u,v;
 double mu,a,b,de,uu,vv,wu,wv,uv; 
 static vector dd[5]={{0.0,0.0,1.0},{0.0,1.0,0.0},{-1.0,0.0,0.0},{0.0,-1.0,0.0},{0.0,0.0,1.0}};
 VECCOPY(dd[index],d)  // direction 
 mu=DOT(d,n);
 if(fabs(mu) < 0.0001)return FALSE;   /* parallel to face */
 VECSUB(pf0,p1,v1)
 mu=DOT(n,v1)/mu;
 if(mu <= 0.0)return FALSE;          /* behind            */
 if(mu >= len)return FALSE;          /* beyond            */
 VECSCALE(mu,d,v1)
 VECSUM(p1,v1,pi)   // intersect point
 VECSUB(pi,pf0,w)
 VECSUB(pf1,pf0,u)
 VECSUB(pf2,pf0,v)
 uv=DOT(u,v);
 uu=DOT(u,u);
 vv=DOT(v,v);
 de=uu*vv-uv*uv;
 if(fabs(de) < 32)return FALSE;  // smallest cell
 wu=DOT(w,u);
 wv=DOT(w,v);
 a=(wu*vv-wv*uv)/de; 
 b=(wv*uu-wu*uv)/de;
 if( a <= 0.0 || a >= 1.0 || b <= 0.0|| b >= 1.0 || (a+b) >= 1.0)return FALSE;
 return TRUE; // strictly inside
} 

static BOOL TouchingFaceOutsideBox(node *n, ivector ip1, ivector ip2, ivector ip3){
 long ma,i;
 for(i=0;i<3;i++){ 
   if(ip1[i] == n->ma[i] && ip2[i] >= n->ma[i] && ip3[i] >= n->ma[i])return TRUE;
   if(ip1[i] == n->mi[i] && ip2[i] <= n->mi[i] && ip3[i] <= n->mi[i])return TRUE;
 }
 return FALSE;
}

static BOOL FaceInVoxel(node *n, tvertex *vBase, tface *f){
 // does face "f" lie inside node "n" 
 long i,j,c,c0,c1,c2,ma,mi;
 static long vxi1[12]={0,1,2,3,4,5,6,7,0,1,2,3}, // link from this voxel vertex 
             vxi2[12]={1,2,3,0,5,6,7,4,4,5,6,7}, // to this voxel vertex 
             vxdi[12]={0,1,2,3,0,1,2,3,4,4,4,4}, // in this directions (0=x 1=y= 2=z)
             len[12];
 tvertex *v;  
 ivector p1,p2,voxv[8];
 // it will return false if every vertex (3) is on the outside of each of the (6) faces
 // eg. if every vertex x coord is > max x coord
 for(i=0;i<3;i++){ // coords
   ma=n->ma[i]; mi=n->mi[i];
   c0=(vBase+f->V[0])->ip[i]; c1=(vBase+f->V[1])->ip[i]; c2=(vBase+f->V[2])->ip[i];
   if(c0 > ma && c1 > ma && c2 > ma)return FALSE;
   if(c0 < mi && c1 < mi && c2 < mi)return FALSE;
 }
 // it will return true if any vertex is strictly inside (3)
 for(j=0;j<3;j++){ 
   v=(vBase+f->V[j]);
   if(v->ip[0] > n->mi[0] && v->ip[0] < n->ma[0] &&
      v->ip[1] > n->mi[1] && v->ip[1] < n->ma[1] &&
      v->ip[2] > n->mi[2] && v->ip[2] < n->ma[2])return TRUE;
 } 
 // it will return true if every vertex (3)  is on a voxel face (6)
 for(i=0;i<6;i++){  
   if(OnFaceEdge(n,i,(vBase+f->V[0])->ip) && 
      OnFaceEdge(n,i,(vBase+f->V[1])->ip) && 
      OnFaceEdge(n,i,(vBase+f->V[2])->ip)) return TRUE;
 }
 // if will return false if one vertex is on a voxel face edge plane and the other two are
 // on the ourside of that face plane
 if(TouchingFaceOutsideBox(n,(vBase+f->V[0])->ip,(vBase+f->V[1])->ip,(vBase+f->V[2])->ip))return FALSE;
 if(TouchingFaceOutsideBox(n,(vBase+f->V[1])->ip,(vBase+f->V[2])->ip,(vBase+f->V[0])->ip))return FALSE;
 if(TouchingFaceOutsideBox(n,(vBase+f->V[2])->ip,(vBase+f->V[0])->ip,(vBase+f->V[1])->ip))return FALSE;
 // it will return true if any face edge (3) strictly crossed one of the (6) voxel faces
 for(i=0;i<3;i++){
   if     (i == 0){ VECCOPY((vBase+f->V[0])->ip,p1)  VECCOPY((vBase+f->V[1])->ip,p2) }
   else if(i == 1){ VECCOPY((vBase+f->V[1])->ip,p1)  VECCOPY((vBase+f->V[2])->ip,p2) }
   else           { VECCOPY((vBase+f->V[2])->ip,p1)  VECCOPY((vBase+f->V[0])->ip,p2) }
   if(EdgeCrossBoxFaceX(p1,p2,n->mi[0],  n->mi[1],n->ma[1],n->mi[2],n->ma[2]))return TRUE;
   if(EdgeCrossBoxFaceX(p1,p2,n->ma[0],  n->mi[1],n->ma[1],n->mi[2],n->ma[2]))return TRUE;
   if(EdgeCrossBoxFaceY(p1,p2,n->mi[1],  n->mi[2],n->ma[2],n->mi[0],n->ma[0]))return TRUE;
   if(EdgeCrossBoxFaceY(p1,p2,n->ma[1],  n->mi[2],n->ma[2],n->mi[0],n->ma[0]))return TRUE;
   if(EdgeCrossBoxFaceZ(p1,p2,n->mi[2],  n->mi[0],n->ma[0],n->mi[1],n->ma[1]))return TRUE;
   if(EdgeCrossBoxFaceZ(p1,p2,n->ma[2],  n->mi[0],n->ma[0],n->mi[1],n->ma[1]))return TRUE;
 }
 // it will return true if any voxel edge (12) passes through a face
 voxv[0][0]=n->mi[0]; voxv[0][1]=n->mi[1]; voxv[0][2]=n->mi[2];
 voxv[1][0]=n->ma[0]; voxv[1][1]=n->mi[1]; voxv[1][2]=n->mi[2];
 voxv[2][0]=n->ma[0]; voxv[2][1]=n->ma[1]; voxv[2][2]=n->mi[2];
 voxv[3][0]=n->mi[0]; voxv[3][1]=n->ma[1]; voxv[3][2]=n->mi[2];
 voxv[4][0]=n->mi[0]; voxv[4][1]=n->mi[1]; voxv[4][2]=n->ma[2];
 voxv[5][0]=n->ma[0]; voxv[5][1]=n->mi[1]; voxv[5][2]=n->ma[2];
 voxv[6][0]=n->ma[0]; voxv[6][1]=n->ma[1]; voxv[6][2]=n->ma[2];
 voxv[7][0]=n->mi[0]; voxv[7][1]=n->ma[1]; voxv[7][2]=n->ma[2];
 len[0]=len[2]=len[4]=len[6]=n->ma[0]-n->mi[0];
 len[1]=len[3]=len[5]=len[7]=n->ma[1]-n->mi[1];
 len[8]=len[9]=len[10]=len[11]=n->ma[2]-n->mi[2];
 for(i=0;i<12;i++){
   if(BoxEdgeCrossFace(voxv[vxi1[i]],voxv[vxi2[i]],vxdi[i],len[i],f->n,
      (vBase+f->V[0])->ip,(vBase+f->V[1])->ip,(vBase+f->V[2])->ip))return TRUE;
 }
 return FALSE;
}

static BOOL PartitionVoxel(node *np, node *na, node *nb, long p, long index){ // index = 2 for z coord 0 for x coord
  // split polygons in node "np" into the two sub nodes "na" and "nb" at plane "p" 
 long i,nnp,nx,n1,n2,o,f,c,sz,xvf,nl,nr,ns,ng;
 tface *F;
 tvertex *V,*v0,*v1,*v2;
 nnp=np->n;
 // assign max and min for the new voxels
 //if(debug != NULL)fprintf(debug,"Split voxel with %ld faces in plane %ld at %ld\n",nnp,index,p); 
 for(i=0;i<3;i++){
   na->mi[i]=nb->mi[i]=np->mi[i];
   na->ma[i]=nb->ma[i]=np->ma[i];
 }
 na->ma[index]=p;     //  "na" points at voxel with the lesser coord value  
 nb->mi[index]=p;     //  "nb" to one with greater coord value
 // calculate number of extra vertices and faces required  
 for(i=0;i<ObjectCount;i++)ObjectXlist[i]=0;
 for(i=0;i<nnp;i++){  
   *(np->vxi+i) = 0;
 }
 // allocate space for extra faces and vertices - this is absolute
 for(i=0,nx=0;i<ObjectCount;i++)if((xvf=ObjectXlist[i]) > 0){
   sz=(Object[i].NoTraceVertices+xvf)*(long)sizeof(tvertex);
   if((V=(tvertex  *)X__Realloc(Object[i].Vtracebase,sz)) == NULL)goto FAILED;
   Object[i].Vtracebase=V;
   sz=(Object[i].NoTraceFaces+xvf)*(long)sizeof(tface);
   if((F=(tface  *)X__Realloc(Object[i].Ftracebase,sz)) == NULL)goto FAILED;
   Object[i].Ftracebase=F;
   nx += xvf;
 }
 // allocate space for list of vertices and in child nodes. This is relative to voxel.
 // this will be over estimate since some of the added faces will go in one voxel and
 // some in the other
 ng=nnp+nx;
 //fprintf(debug,"Guessed voxel count %ld\n",ng);
 if(ng > 0){
   if((na->oid=(UC *)X__Malloc(sizeof(UC)*ng)) == NULL)goto FAILED;
   if((na->fid=(US *)X__Malloc(sizeof(US)*ng)) == NULL)goto FAILED;
   if((na->vxi=(UC *)X__Malloc(sizeof(UC)*ng)) == NULL)goto FAILED;
   if((nb->oid=(UC *)X__Malloc(sizeof(UC)*ng)) == NULL)goto FAILED;
   if((nb->fid=(US *)X__Malloc(sizeof(US)*ng)) == NULL)goto FAILED;
   if((nb->vxi=(UC *)X__Malloc(sizeof(UC)*ng)) == NULL)goto FAILED;
   for(i=0;i<ng;i++){na->oid[i]=999,na->fid[i]=999,na->vxi[i]=999;nb->oid[i]=666,nb->fid[i]=666,nb->vxi[i]=666;} // For testing
 }
 // split vertices and faces - allocate them to the appropriate sub voxel - use VXI code
 nl=nr=ns=0;
 for(i=0,n1=0,n2=0;i<nnp;i++){  
   c = (long)(*(np->vxi+i));   
   o = (long)(*(np->oid+i));
   f = (long)(*(np->fid+i));
//index 0 = x plane
   F=(Object[o].Ftracebase+f);
   V=Object[o].Vtracebase;
   if(FaceInVoxel(na,V,F)){
     *(na->fid+n1)=f; 
     *(na->oid+n1)=o; n1++;
     nl++;
   }
   if(FaceInVoxel(nb,V,F)) {
     *(nb->fid+n2)=f; 
     *(nb->oid+n2)=o; n2++;
     nr++;
   }
 }
 //fprintf(debug,"Assigned to lower %ld  to upper %ld  and split %ld\n",nl,nr,ns);
 // assign the correct sizes to voxel face list
 if(ng > 0 && n1 == 0){
   X__Free(na->oid); na->oid=NULL;
   X__Free(na->fid); na->fid=NULL;
   X__Free(na->vxi); na->vxi=NULL;
 }
 else if(n1 > 0){
   na->oid=(UC *)X__Realloc(na->oid,sizeof(UC)*n1);
   na->vxi=(UC *)X__Realloc(na->vxi,sizeof(UC)*n1);
   na->fid=(US *)X__Realloc(na->fid,sizeof(US)*n1);
 }
 if(ng > 0 && n2 == 0){
   X__Free(nb->oid); nb->oid=NULL;
   X__Free(nb->fid); nb->fid=NULL;
   X__Free(nb->vxi); nb->vxi=NULL;
 }
 if(n2 > 0){
   nb->oid=(UC *)X__Realloc(nb->oid,sizeof(UC)*n2);
   nb->vxi=(UC *)X__Realloc(nb->vxi,sizeof(UC)*n2);
   nb->fid=(US *)X__Realloc(nb->fid,sizeof(US)*n2);
 }
 // parent voxel now split so free face list
 if(np->oid != NULL)X__Free(np->oid); np->oid=NULL;
 if(np->fid != NULL)X__Free(np->fid); np->fid=NULL;
 if(np->vxi != NULL)X__Free(np->vxi); np->vxi=NULL;
 na->n=n1;
 nb->n=n2;
 //fprintf(debug,"%ld split into %ld  %ld\n",np->n, na->n,nb->n); fflush(debug); 
 np->n=0; 
 return OK;
 FAILED:
 return FAIL;
}