Logo Search packages:      
Sourcecode: pymol version File versions  Download package

Ray.c

/* 
A* -------------------------------------------------------------------
B* This file contains source code for the PyMOL computer program
C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific. 
D* -------------------------------------------------------------------
E* It is unlawful to modify or remove this copyright notice.
F* -------------------------------------------------------------------
G* Please see the accompanying LICENSE file for further information. 
H* -------------------------------------------------------------------
I* Additional authors of this source file include:
-*   Larry Coopet (various optimizations)
-*   Chris Want (RayRenderVRML2, via the public domain )
-*
Z* -------------------------------------------------------------------
*/


#include"os_predef.h"
#include"os_std.h"

#include"Base.h"
#include"MemoryDebug.h"
#include"Err.h"
#include"Vector.h"
#include"OOMac.h"
#include"Setting.h"
#include"Ortho.h"
#include"Util.h"
#include"Ray.h"
#include"Triangle.h" 
#include"Color.h"
#include"Matrix.h"
#include"P.h"
#include"MemoryCache.h"
#include"Character.h"
#include"Text.h"
#include"PyMOL.h"
#include"Scene.h"
#include"PConv.h"

#ifdef _PYMOL_INLINE
#undef _PYMOL_INLINE
#include"Basis.c"
#define _PYMOL_INLINE
#endif

#ifndef RAY_SMALL
#define RAY_SMALL 0.00001
#endif

/* BASES 
   0 contains untransformed vertices (vector size = 3)
      1 contains transformed vertices (vector size = 3)
      2 contains transformed vertices for shadowing 
*/

/* note: the following value must be at least one greater than the max
   number of lights */

#define MAX_BASIS 12

typedef float float3[3];
typedef float float4[4];

struct _CRayThreadInfo {
  CRay *ray;
  int width,height;
  unsigned int *image;
  float front,back;
  unsigned int fore_mask;
  float *bkrd;
  float ambient;
  unsigned int background;
  int border;
  int phase, n_thread;
  int x_start,x_stop;
  int y_start,y_stop;
  unsigned int *edging;
  unsigned int edging_cutoff;
  int perspective;
  float fov,pos[3];
  float *depth;
};

 struct _CRayHashThreadInfo {
   CBasis *basis;
   int *vert2prim;
   CPrimitive *prim;
   int n_prim;
   float *clipBox;
   unsigned int *image;
   unsigned int background;
   unsigned int bytes;
   int perspective;
   float front;
   int phase;
   float size_hint;
   CRay *ray;
 };

 struct _CRayAntiThreadInfo {
  unsigned int *image;
  unsigned int *image_copy;
  unsigned int width,height;
  int mag;
  int phase,n_thread;
  CRay *ray;
};

void RayRelease(CRay *I);
void RayWobble(CRay *I,int mode,float *v);
void RayTransparentf(CRay *I,float v);

void RaySetup(CRay *I);
void RayColor3fv(CRay *I,float *v);
void RaySphere3fv(CRay *I,float *v,float r);
void RayCharacter(CRay *I,int char_id);
void RayCylinder3fv(CRay *I,float *v1,float *v2,float r,float *c1,float *c2);
void RaySausage3fv(CRay *I,float *v1,float *v2,float r,float *c1,float *c2);
void RayInteriorColor3fv(CRay *I,float *v,int passive);

void RayTriangle3fv(CRay *I,
                                      float *v1,float *v2,float *v3,
                                      float *n1,float *n2,float *n3,
                                      float *c1,float *c2,float *c3);

void RayTriangleTrans3fv(CRay *I,
                         float *v1,float *v2,float *v3,
                         float *n1,float *n2,float *n3,
                         float *c1,float *c2,float *c3,
                         float  t1,float  t2,float  t3);
void RayEllipsoid3fv(CRay *I,
                     float *v,float r,
                     float *n1,float *n2,float *n3);

void RayApplyMatrix33( unsigned int n, float3 *q, const float m[16],
                                          float3 *p );
void RayApplyMatrixInverse33( unsigned int n, float3 *q, const float m[16],
                              float3 *p );

void RayExpandPrimitives(CRay *I);
void RayTransformBasis(CRay *I,CBasis *B,int group_id);

int PrimitiveSphereHit(CRay *I,float *v,float *n,float *minDist,int except);

void RayTransformNormals33( unsigned int n, float3 *q, const float m[16],float3 *p );
void RayTransformInverseNormals33( unsigned int n, float3 *q, const float m[16],float3 *p );
void RayProjectTriangle(CRay *I,RayInfo *r,float *light,float *v0,float *n0,float scale);
void RayCustomCylinder3fv(CRay *I,float *v1,float *v2,float r,
                          float *c1,float *c2,int cap1,int cap2);
void RaySetContext(CRay *I,int context)
{
  if(context>=0)
    I->Context=context;
  else
    I->Context=0;
}
void RayApplyContextToNormal(CRay *I,float *v);
void RayApplyContextToVertex(CRay *I,float *v);

static float RayGetScreenVertexScale(CRay *I,float *v1)
{
  /* what size should a screen pixel be at the coordinate provided? */

  float vt[3];
  float front_size,ratio;
  RayApplyMatrix33(1,(float3*)vt,I->ModelView,(float3*)v1);

  front_size = 2*I->Volume[4]*((float)tan((I->Fov/2.0F)*PI/180.0F))/(I->Height);
  ratio = front_size*(-vt[2]/I->Volume[4]);
  return ratio;
}

void RayApplyContextToVertex(CRay *I,float *v)
{
  switch(I->Context) {
  case 1:
    {
      float tw;
      float th;
      
      if(I->AspRatio>1.0F) {
        tw = I->AspRatio;
        th = 1.0F;
      } else {
        th = 1.0F/I->AspRatio;
        tw = 1.0F;
      }
      if(!SettingGetGlobal_b(I->G,cSetting_ortho)) {
        float scale = v[2]+0.5F;
        scale = I->FrontBackRatio*scale + 1.0F - scale;
          
        /* z-coodinate is easy... */

        v[2]=v[2]*I->Range[2]-(I->Volume[4]+I->Volume[5])/2.0F;
        v[0]-=0.5F;
        v[1]-=0.5F;
        v[0]=scale*v[0]*I->Range[0]/tw+(I->Volume[0]+I->Volume[1])/2.0F; 
        v[1]=scale*v[1]*I->Range[1]/th+(I->Volume[2]+I->Volume[3])/2.0F; 
        
        RayApplyMatrixInverse33(1,(float3*)v,I->ModelView,(float3*)v);    
      } else {
        v[0]+=(tw-1.0F)/2;
        v[1]+=(th-1.0F)/2;
        v[0]=v[0]*(I->Range[0]/tw)+I->Volume[0];
        v[1]=v[1]*(I->Range[1]/th)+I->Volume[2];
        v[2]=v[2]*I->Range[2]-(I->Volume[4]+I->Volume[5])/2.0F;
        RayApplyMatrixInverse33(1,(float3*)v,I->ModelView,(float3*)v);    
      }


    }
    break;
  }
}
void RayApplyContextToNormal(CRay *I,float *v)
{
  switch(I->Context) {
  case 1:
    RayTransformInverseNormals33(1,(float3*)v,I->ModelView,(float3*)v);    
    break;
  }
}
int RayGetNPrimitives(CRay *I)    
{
  return(I->NPrimitive);
}
/*========================================================================*/

#ifdef _PYMOL_INLINE
__inline__
#endif
static void RayGetSphereNormal(CRay *I,RayInfo *r)
{
  
  r->impact[0]=r->base[0]; 
  r->impact[1]=r->base[1]; 
  r->impact[2]=r->base[2]-r->dist;
  
  r->surfnormal[0]=r->impact[0]-r->sphere[0];
  r->surfnormal[1]=r->impact[1]-r->sphere[1];
  r->surfnormal[2]=r->impact[2]-r->sphere[2];
  
  normalize3f(r->surfnormal);
}

#ifdef _PYMOL_INLINE
__inline__
#endif
static void RayGetSphereNormalPerspective(CRay *I,RayInfo *r)
{
  
  r->impact[0]=r->base[0] + r->dist*r->dir[0];
  r->impact[1]=r->base[1] + r->dist*r->dir[1];
  r->impact[2]=r->base[2] + r->dist*r->dir[2];
  
  r->surfnormal[0]=r->impact[0]-r->sphere[0];
  r->surfnormal[1]=r->impact[1]-r->sphere[1];
  r->surfnormal[2]=r->impact[2]-r->sphere[2];
  
  normalize3f(r->surfnormal);
}

static void fill(unsigned int *buffer, unsigned int value,unsigned int cnt)
{
  while(cnt&0xFFFFFF80) {
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    cnt-=0x20;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
    *(buffer++) = value;
  }
  while(cnt--) {
    *(buffer++) = value;
  }
}
/*========================================================================*/
#ifdef _PYMOL_INLINE
__inline__
#endif
static void RayReflectAndTexture(CRay *I,RayInfo *r,int perspective)
{

  if(r->prim->wobble)
    switch(r->prim->wobble) {
    case 1:
      scatter3f(r->surfnormal,I->WobbleParam[0]);
      break;
    case 2:
      wiggle3f(r->surfnormal,r->impact,I->WobbleParam);
      break;
    case 3: 
      {
        float3 v;
        float3 n;
        copy3f(r->impact,v);
        RayApplyMatrixInverse33(1,&v,I->ModelView,&v);
        n[0]=(float)cos((v[0]+v[1]+v[2])*I->WobbleParam[1]);
        n[1]=(float)cos((v[0]-v[1]+v[2])*I->WobbleParam[1]);
        n[2]=(float)cos((v[0]+v[1]-v[2])*I->WobbleParam[1]);
        RayTransformNormals33(1,&n,I->ModelView,&n);
        scale3f(n,I->WobbleParam[0],n);
        add3f(n,r->surfnormal,r->surfnormal);
        normalize3f(r->surfnormal);
      }
    case 4: 
      {
        float3 v;
        float3 n;
        float *tp = I->WobbleParam;
        copy3f(r->impact,v);
        RayApplyMatrixInverse33(1,&v,I->ModelView,&v);
        n[0]=I->Random[0xFF&(int)((cos((v[0])*tp[1])*256*tp[2]))];
        n[1]=I->Random[0xFF&(int)((cos((v[1])*tp[1])*256*tp[2]+96))];
        n[2]=I->Random[0xFF&(int)((cos((v[2])*tp[1])*256*tp[2]+148))];
        RayTransformNormals33(1,&n,I->ModelView,&n);
        scale3f(n,tp[0],n);
        add3f(n,r->surfnormal,r->surfnormal);
        normalize3f(r->surfnormal);
      }
      break;
    case 5: 
      {
        float3 v;
        float3 n;
        float *tp = I->WobbleParam;
        copy3f(r->impact,v);
        RayApplyMatrixInverse33(1,&v,I->ModelView,&v);
        n[0]=I->Random[0xFF&(int)((v[0]*tp[1])+0)]+
          I->Random[0xFF&(int)((v[1]*tp[1])+20)]+
          I->Random[0xFF&(int)((v[2]*tp[1])+40)];
        n[1]=I->Random[0xFF&(int)((-v[0]*tp[1])+90)]+
          I->Random[0xFF&(int)((v[1]*tp[1])+100)]+
          I->Random[0xFF&(int)((-v[2]*tp[1])+120)];
        n[2]=I->Random[0xFF&(int)((v[0]*tp[1])+200)]+
          I->Random[0xFF&(int)((-v[1]*tp[1])+70)]+
          I->Random[0xFF&(int)((v[2]*tp[1])+30)];
        
        n[0]+=
          I->Random[0xFF&((int)((v[0]-v[1])*tp[1])+0)] +
          I->Random[0xFF&((int)((v[1]-v[2])*tp[1])+20)] +
          I->Random[0xFF&((int)((v[2]-v[0])*tp[1])+40)];
        n[1]+=
          I->Random[0xFF&((int)((v[0]+v[1])*tp[1])+10)]+
          I->Random[0xFF&((int)((v[1]+v[2])*tp[1])+90)]+
          I->Random[0xFF&((int)((v[2]+v[0])*tp[1])+30)];
        n[2]+=
          I->Random[0xFF&((int)((-v[0]+v[1])*tp[1])+220)]+
          I->Random[0xFF&((int)((-v[1]+v[2])*tp[1])+20)]+
          I->Random[0xFF&((int)((-v[2]+v[0])*tp[1])+50)];
        
        n[0]+=
          I->Random[0xFF&((int)((v[0]+v[1]+v[2])*tp[1])+5)]+
          I->Random[0xFF&((int)((v[0]+v[1]+v[2])*tp[1])+25)]+
          I->Random[0xFF&((int)((v[0]+v[1]+v[2])*tp[1])+46)];
        n[1]+=
          I->Random[0xFF&((int)((-v[0]-v[1]+v[2])*tp[1])+90)]+
          I->Random[0xFF&((int)((-v[0]-v[1]+v[2])*tp[1])+45)]+
          I->Random[0xFF&((int)((-v[0]-v[1]+v[2])*tp[1])+176)];
        n[2]+=
          I->Random[0xFF&((int)((v[0]+v[1]-v[2])*tp[1])+192)]+
          I->Random[0xFF&((int)((v[0]+v[1]-v[2])*tp[1])+223)]+
          I->Random[0xFF&((int)((v[0]+v[1]-v[2])*tp[1])+250)];

        RayTransformNormals33(1,&n,I->ModelView,&n);
        scale3f(n,tp[0],n);
        add3f(n,r->surfnormal,r->surfnormal);
        normalize3f(r->surfnormal);
      }
      break;
    }
  if(perspective) {
    r->dotgle = dot_product3f(r->dir, r->surfnormal);
    r->flat_dotgle = -r->dotgle;
  
    r->reflect[0]= r->dir[0] - ( 2 * r->dotgle * r->surfnormal[0] );
    r->reflect[1]= r->dir[1] - ( 2 * r->dotgle * r->surfnormal[1] );
    r->reflect[2]= r->dir[2] - ( 2 * r->dotgle * r->surfnormal[2] );
  } else {
    r->dotgle = -r->surfnormal[2]; 
    r->flat_dotgle = r->surfnormal[2];
    
    r->reflect[0]= - ( 2 * r->dotgle * r->surfnormal[0] );
    r->reflect[1]= - ( 2 * r->dotgle * r->surfnormal[1] );
    r->reflect[2]= -1.0F - ( 2 * r->dotgle * r->surfnormal[2] );
  }
}
/*========================================================================*/
void RayExpandPrimitives(CRay *I)
{
  int a;
  float *v0,*v1,*n0,*n1;
  CBasis *basis;
  int nVert, nNorm;
  float voxel_floor;

  nVert=0;
  nNorm=0;
  for(a=0;a<I->NPrimitive;a++) {
       switch(I->Primitive[a].type) {
       case cPrimSphere:
       nVert++;
       break;
     case cPrimEllipsoid:
       nVert++;
       nNorm+=3;
       break;
       case cPrimCylinder:
     case cPrimSausage:
       nVert++;
       nNorm++;
       break;
       case cPrimTriangle:
       case cPrimCharacter:
       nVert+=3;
       nNorm+=4;
       break;
       }
  }

  basis = I->Basis;
  
  VLACacheSize(I->G,basis->Vertex,float,3*nVert,0,cCache_basis_vertex);
  VLACacheSize(I->G,basis->Radius,float,nVert,0,cCache_basis_radius);
  VLACacheSize(I->G,basis->Radius2,float,nVert,0,cCache_basis_radius2);
  VLACacheSize(I->G,basis->Vert2Normal,int,nVert,0,cCache_basis_vert2normal);
  VLACacheSize(I->G,basis->Normal,float,3*nNorm,0,cCache_basis_normal);

  VLACacheSize(I->G,I->Vert2Prim,int,nVert,0,cCache_ray_vert2prim);

  voxel_floor=I->PixelRadius/2.0F;

  basis->MaxRadius = 0.0F;
  basis->MinVoxel = 0.0F;
  basis->NVertex=nVert;
  basis->NNormal=nNorm;

  nVert=0;
  nNorm=0;
  v0=basis->Vertex;
  n0=basis->Normal;
  for(a=0;a<I->NPrimitive;a++) {
       switch(I->Primitive[a].type) {
       case cPrimTriangle:
       case cPrimCharacter:

            I->Primitive[a].vert=nVert;
            I->Vert2Prim[nVert]=a;
            I->Vert2Prim[nVert+1]=a;
            I->Vert2Prim[nVert+2]=a;
            basis->Radius[nVert]=I->Primitive[a].r1;
            basis->Radius2[nVert]=I->Primitive[a].r1*I->Primitive[a].r1; /*necessary??*/
            /*          if(basis->Radius[nVert]>basis->MinVoxel)
                        basis->MinVoxel=basis->Radius[nVert];*/
            if(basis->MinVoxel<voxel_floor)
              basis->MinVoxel=voxel_floor;
            basis->Vert2Normal[nVert]=nNorm;
            basis->Vert2Normal[nVert+1]=nNorm;
            basis->Vert2Normal[nVert+2]=nNorm;
            n1=I->Primitive[a].n0;
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            n1=I->Primitive[a].n1;
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            n1=I->Primitive[a].n2;
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            n1=I->Primitive[a].n3;
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            nNorm+=4;
            v1=I->Primitive[a].v1;
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            v1=I->Primitive[a].v2;
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            v1=I->Primitive[a].v3;
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            nVert+=3;
            break;
       case cPrimSphere:
            I->Primitive[a].vert=nVert;
            I->Vert2Prim[nVert]=a;
            v1=I->Primitive[a].v1;
            basis->Radius[nVert]=I->Primitive[a].r1;
            basis->Radius2[nVert]=I->Primitive[a].r1*I->Primitive[a].r1; /*precompute*/
            if(basis->Radius[nVert]>basis->MaxRadius)
              basis->MaxRadius=basis->Radius[nVert];
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            nVert++;
        break;
       case cPrimEllipsoid:
            I->Primitive[a].vert=nVert;
            I->Vert2Prim[nVert]=a;
            v1=I->Primitive[a].v1;
            basis->Radius[nVert]=I->Primitive[a].r1;
            basis->Radius2[nVert]=I->Primitive[a].r1*I->Primitive[a].r1; /*precompute*/
            if(basis->Radius[nVert]>basis->MaxRadius)
              basis->MaxRadius=basis->Radius[nVert];
            basis->Vert2Normal[nVert]=nNorm;
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            nVert++;
            n1=I->Primitive[a].n1;
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            n1=I->Primitive[a].n2;
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            n1=I->Primitive[a].n3;
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            (*n0++)=(*n1++);
            nNorm+=3;
            break;
       case cPrimCylinder:
       case cPrimSausage:
            I->Primitive[a].vert=nVert;
            I->Vert2Prim[nVert]=a;
            basis->Radius[nVert]=I->Primitive[a].r1;
            basis->Radius2[nVert]=I->Primitive[a].r1*I->Primitive[a].r1; /*precompute*/
            if(basis->MinVoxel<voxel_floor)
        basis->MinVoxel=voxel_floor;
            subtract3f(I->Primitive[a].v2,I->Primitive[a].v1,n0);
            I->Primitive[a].l1=(float)length3f(n0);
            normalize3f(n0);
            n0+=3;
            basis->Vert2Normal[nVert]=nNorm;
            nNorm++;
            v1=I->Primitive[a].v1;
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            (*v0++)=(*v1++);
            nVert++;
            break;
       }
  }
  if(nVert>basis->NVertex) {
    fprintf(stderr,"Error: basis->NVertex exceeded\n");
  }
  PRINTFB(I->G,FB_Ray,FB_Blather)
    " Ray: minvoxel  %8.3f\n Ray: NPrimit  %d nvert %d\n",basis->MinVoxel,I->NPrimitive,nVert
    ENDFB(I->G);
}
/*========================================================================*/
static void RayComputeBox(CRay *I)
{

#define minmax(v,r) { \
  xp = v[0] + r;\
  xm = v[0] - r;\
  yp = v[1] + r;\
  ym = v[1] - r;\
  zp = v[2] + r;\
  zm = v[2] - r;\
  if(xmin>xm) xmin = xm;\
  if(xmax<xp) xmax = xp;\
  if(ymin>ym) ymin = ym;\
  if(ymax<yp) ymax = yp;\
  if(zmin>zm) zmin = zm;\
  if(zmax<zp) zmax = zp;\
}

  CPrimitive *prm;
  CBasis *basis1;

  float xmin=0.0F,ymin=0.0F,xmax=0.0F,ymax=0.0F,zmin=0.0F,zmax=0.0F;
  float xp,xm,yp,ym,zp,zm;

  float *v,r;
  float vt[3];
  const float _0 = 0.0F;
  int a;

  basis1 = I->Basis+1;
  if(basis1->NVertex) {
    xmin = xmax = basis1->Vertex[0];
    ymin = ymax = basis1->Vertex[1];
    zmin = zmax = basis1->Vertex[2];

    for(a=0;a<I->NPrimitive;a++) {
      prm=I->Primitive+a;
      
      switch(prm->type) 
      {
      case cPrimTriangle:
      case cPrimCharacter:

        r = _0;
        v = basis1->Vertex + prm->vert*3;
        minmax(v,r);
        v = basis1->Vertex + prm->vert*3+3;
        minmax(v,r);
        v = basis1->Vertex + prm->vert*3+6;
        minmax(v,r);
        break;
      case cPrimSphere:
      case cPrimEllipsoid:
        r = prm->r1;
        v = basis1->Vertex + prm->vert*3;
        minmax(v,r);
        break;
      case cPrimCylinder:
      case cPrimSausage:
        r = prm->r1;
        v = basis1->Vertex + prm->vert*3;
        minmax(v,r);
        v = basis1->Normal+basis1->Vert2Normal[prm->vert]*3;
        scale3f(v,prm->l1,vt);
        v = basis1->Vertex + prm->vert*3;
        add3f(v,vt,vt);
        minmax(vt,r);
        break;
      }     /* end of switch */
       }
  }
  I->min_box[0] = xmin;
  I->min_box[1] = ymin;
  I->min_box[2] = zmin;
  I->max_box[0] = xmax;
  I->max_box[1] = ymax;
  I->max_box[2] = zmax;
}

static void RayTransformFirst(CRay *I,int perspective)
{
  CBasis *basis0,*basis1;
  CPrimitive *prm;
  int a;
  float *v0;
  int backface_cull;

  backface_cull = (int)SettingGet(I->G,cSetting_backface_cull);
  
  if((SettingGet(I->G,cSetting_two_sided_lighting)||
      (SettingGet(I->G,cSetting_transparency_mode)==1)||
      (SettingGet(I->G,cSetting_ray_interior_color)!=-1)||
      I->CheckInterior))
     backface_cull=0;

  basis0 = I->Basis;
  basis1 = I->Basis+1;
  
  VLACacheSize(I->G,basis1->Vertex,float,3*basis0->NVertex,1,cCache_basis_vertex);
  VLACacheSize(I->G,basis1->Normal,float,3*basis0->NNormal,1,cCache_basis_normal);
  VLACacheSize(I->G,basis1->Precomp,float,3*basis0->NNormal,1,cCache_basis_precomp);
  VLACacheSize(I->G,basis1->Vert2Normal,int,basis0->NVertex,1,cCache_basis_vert2normal);
  VLACacheSize(I->G,basis1->Radius,float,basis0->NVertex,1,cCache_basis_radius);
  VLACacheSize(I->G,basis1->Radius2,float,basis0->NVertex,1,cCache_basis_radius2);
  
  RayApplyMatrix33(basis0->NVertex,(float3*)basis1->Vertex,
                                I->ModelView,(float3*)basis0->Vertex);

  for(a=0;a<basis0->NVertex;a++) {
    basis1->Radius[a]=basis0->Radius[a];
    basis1->Radius2[a]=basis0->Radius2[a];
    basis1->Vert2Normal[a]=basis0->Vert2Normal[a];
  }
  basis1->MaxRadius=basis0->MaxRadius;
  basis1->MinVoxel=basis0->MinVoxel;
  basis1->NVertex=basis0->NVertex;

  RayTransformNormals33(basis0->NNormal,(float3*)basis1->Normal,
                                I->ModelView,(float3*)basis0->Normal);
  
  basis1->NNormal=basis0->NNormal;

  if(perspective) {
    for(a=0;a<I->NPrimitive;a++) {
      prm=I->Primitive+a;
      
      prm=I->Primitive+a;
      switch(prm->type) {
      case cPrimTriangle:
      case cPrimCharacter:
        BasisTrianglePrecomputePerspective(
                                           basis1->Vertex+prm->vert*3,
                                           basis1->Vertex+prm->vert*3+3,
                                           basis1->Vertex+prm->vert*3+6,
                                           basis1->Precomp+basis1->Vert2Normal[prm->vert]*3);
        break;
      }
    }
  } else {
    for(a=0;a<I->NPrimitive;a++) {
      prm=I->Primitive+a;
      switch(prm->type) {
      case cPrimTriangle:
      case cPrimCharacter:
        
        BasisTrianglePrecompute(basis1->Vertex+prm->vert*3,
                                basis1->Vertex+prm->vert*3+3,
                                basis1->Vertex+prm->vert*3+6,
                                basis1->Precomp+basis1->Vert2Normal[prm->vert]*3);
        v0 = basis1->Normal + (basis1->Vert2Normal[prm->vert]*3 + 3);
        prm->cull = backface_cull&&((v0[2]<0.0F)&&(v0[5]<0.0F)&&(v0[8]<0.0F));
        break;
      case cPrimSausage:
      case cPrimCylinder:
        BasisCylinderSausagePrecompute(basis1->Normal + basis1->Vert2Normal[prm->vert]*3,
                                       basis1->Precomp + basis1->Vert2Normal[prm->vert]*3);
        break;
        
      }
    }
  }
}
/*========================================================================*/
void RayTransformBasis(CRay *I,CBasis *basis1,int group_id)
{
  CBasis *basis0;
  int a;
  float *v0,*v1;
  CPrimitive *prm;

  basis0 = I->Basis+1;

  VLACacheSize(I->G,basis1->Vertex,float,3*basis0->NVertex,group_id,cCache_basis_vertex);
  VLACacheSize(I->G,basis1->Normal,float,3*basis0->NNormal,group_id,cCache_basis_normal);
  VLACacheSize(I->G,basis1->Precomp,float,3*basis0->NNormal,group_id,cCache_basis_precomp);
  VLACacheSize(I->G,basis1->Vert2Normal,int,basis0->NVertex,group_id,cCache_basis_vert2normal);
  VLACacheSize(I->G,basis1->Radius,float,basis0->NVertex,group_id,cCache_basis_radius);
  VLACacheSize(I->G,basis1->Radius2,float,basis0->NVertex,group_id,cCache_basis_radius2);
  v0=basis0->Vertex;
  v1=basis1->Vertex;
  for(a=0;a<basis0->NVertex;a++) {
    matrix_transform33f3f(basis1->Matrix,v0,v1);
    v0+=3;
    v1+=3;
    basis1->Radius[a]=basis0->Radius[a];
    basis1->Radius2[a]=basis0->Radius2[a];
    basis1->Vert2Normal[a]=basis0->Vert2Normal[a];
  }
  v0=basis0->Normal;
  v1=basis1->Normal;
  for(a=0;a<basis0->NNormal;a++) {
    matrix_transform33f3f(basis1->Matrix,v0,v1);
    v0+=3;
    v1+=3;
  }
  basis1->MaxRadius=basis0->MaxRadius;
  basis1->MinVoxel=basis0->MinVoxel;
  basis1->NVertex=basis0->NVertex;
  basis1->NNormal=basis0->NNormal;


  for(a=0;a<I->NPrimitive;a++) {
       prm=I->Primitive+a;
    switch(prm->type) {
    case cPrimTriangle:
       case cPrimCharacter:

        BasisTrianglePrecompute(basis1->Vertex+prm->vert*3,
                                basis1->Vertex+prm->vert*3+3,
                                basis1->Vertex+prm->vert*3+6,
                                basis1->Precomp+basis1->Vert2Normal[prm->vert]*3);
        break;
    case cPrimSausage:
    case cPrimCylinder:
      BasisCylinderSausagePrecompute(basis1->Normal + basis1->Vert2Normal[prm->vert]*3,
                                     basis1->Precomp + basis1->Vert2Normal[prm->vert]*3);
      break;
      
       }
  }
}

/*========================================================================*/
void RayRenderTest(CRay *I,int width,int height,float front,float back,float fov)
{

  PRINTFB(I->G,FB_Ray,FB_Details)
    " RayRenderTest: obtained %i graphics primitives.\n",I->NPrimitive 
    ENDFB(I->G);
}
/*========================================================================*/

G3dPrimitive *RayRenderG3d(CRay *I,int width, int height,
                           float front, float back, float fov,int quiet)
{
  /* generate a rendering stream for Miguel's G3d java rendering engine */

  register float scale_x,scale_y,scale_z;
  int shift_x,shift_y;
  float *d;
  CBasis *base;
  CPrimitive *prim;
  float *vert;
  float vert2[3];
  int a;
  G3dPrimitive *jprim = VLAlloc(G3dPrimitive,10000),*jp;
  int n_jp = 0;

#define convert_r(r) 2*(int)(r*scale_x);
#define convert_x(x) shift_x + (int)(x*scale_x);
#define convert_y(y) height - (shift_y + (int)(y*scale_y));
#define convert_z(z) -(int)((z+front)*scale_x);
#define convert_col(c) (0xFF000000 | (((int)(c[0]*255.0))<<16) | (((int)(c[1]*255.0))<<8) | (((int)(c[2]*255.0))))

  RayExpandPrimitives(I);
  RayTransformFirst(I,0);

  if(!quiet) {
    PRINTFB(I->G,FB_Ray,FB_Details)
      " RayRenderG3d: processed %i graphics primitives.\n",I->NPrimitive 
      ENDFB(I->G);
  }
  base = I->Basis+1;

  /* always orthoscopic */
  
  /* front should give a zero Z, 
     -I->Range[0] should be off the right hand size
     I->Range[1] should be off the top */
  scale_x = width/I->Range[0];
  scale_y = height/I->Range[1];
  scale_z = -4096.0F/(back-front);
  shift_x = width/2;
  shift_y = height/2;

  for(a=0;a<I->NPrimitive;a++) {
    prim = I->Primitive+a;
    vert = base->Vertex+3*(prim->vert);
    switch(prim->type) {
       case cPrimSphere:
      VLACheck(jprim,G3dPrimitive,n_jp);
      jp = jprim + n_jp;
      jp->op = 1;
      jp->r = convert_r(prim->r1);
      jp->x1 = convert_x(vert[0]);
      jp->y1 = convert_y(vert[1]);
      jp->z1 = convert_z(vert[2]);
      jp->c = convert_col(prim->c1);
      n_jp++;
            break;
    case cPrimSausage:
      VLACheck(jprim,G3dPrimitive,n_jp);
      d=base->Normal+3*base->Vert2Normal[prim->vert];
      scale3f(d,prim->l1,vert2);
      add3f(vert,vert2,vert2);

      jp = jprim + n_jp;
      jp->op = 3;
      jp->r = convert_r(prim->r1);
      jp->x1 = convert_x(vert[0]);
      jp->y1 = convert_y(vert[1]);
      jp->z1 = convert_z(vert[2]);
      jp->x2 = convert_x(vert2[0]);
      jp->y2 = convert_y(vert2[1]);
      jp->z2 = convert_z(vert2[2]);
      jp->c = convert_col(prim->c1);
      n_jp++;
            break;
       case cPrimTriangle:
      VLACheck(jprim,G3dPrimitive,n_jp);
      jp = jprim + n_jp;
      jp->op = 2;
      jp->x1 = convert_x(vert[0]);
      jp->y1 = convert_y(vert[1]);
      jp->z1 = convert_z(vert[2]);
      jp->x2 = convert_x(vert[3]);
      jp->y2 = convert_y(vert[4]);
      jp->z2 = convert_z(vert[5]);
      jp->x3 = convert_x(vert[6]);
      jp->y3 = convert_y(vert[7]);
      jp->z3 = convert_z(vert[8]);
      jp->c = convert_col(prim->c1);
      n_jp++;
            break;
    }
  }
  VLASize(jprim,G3dPrimitive,n_jp);
  return jprim;
}

#if 0
void RayRenderVRML1(CRay *I,int width,int height,
                    char **vla_ptr,float front,float back,
                float fov, float angle,float z_corr)
{
  char *vla = *vla_ptr;
  int cc = 0; /* character count */
  OrthoLineType buffer;
  
  RayExpandPrimitives(I);
  RayTransformFirst(I,0);

  strcpy(buffer,"#VRML V1.0 ascii\n\n");
  UtilConcatVLA(&vla,&cc,buffer);

  UtilConcatVLA(&vla,&cc,"MaterialBinding { value OVERALL }\n");

  sprintf(buffer,"Material {\n ambientColor 0 0 0\n diffuseColor 1 1 1\n specularColor 1 1 1\nshininess 0.2\n}\n");
  UtilConcatVLA(&vla,&cc,buffer);

  { 
    int a;
    CPrimitive *prim;
    float *vert;
    CBasis *base = I->Basis+1;

    UtilConcatVLA(&vla,&cc,"Separator {\n");

  UtilConcatVLA(&vla,&cc,"MatrixTransform {\n");
  UtilConcatVLA(&vla,&cc,"matrix 1.0 0.0 0.0 0.0\n");
  UtilConcatVLA(&vla,&cc,"       0.0 1.0 0.0 0.0\n");
  UtilConcatVLA(&vla,&cc,"       0.0 0.0 1.0 0.0\n");
  sprintf(buffer,"    %8.6f %8.6f %8.6f 1.0\n",
          (I->Volume[0]+I->Volume[1])/2,
          (I->Volume[2]+I->Volume[3])/2,0.0F);
  UtilConcatVLA(&vla,&cc,buffer);
  UtilConcatVLA(&vla,&cc,"}\n");

    for(a=0;a<I->NPrimitive;a++) {
      prim = I->Primitive+a;
      vert = base->Vertex+3*(prim->vert);
      switch(prim->type) {
      case cPrimSphere:
        sprintf(buffer,
                "Material {\ndiffuseColor %6.4f %6.4f %6.4f\n}\n\n", 
                prim->c1[0],prim->c1[1],prim->c1[2]);
        UtilConcatVLA(&vla,&cc,buffer);    
        UtilConcatVLA(&vla,&cc,"Separator {\n");
        sprintf(buffer,
                "Transform {\ntranslation %8.6f %8.6f %8.6f\nscaleFactor %8.6f %8.6f %8.6f\n}\n",
                vert[0],vert[1],vert[2]-z_corr,         
                prim->r1,prim->r1,prim->r1);
        UtilConcatVLA(&vla,&cc,buffer);    
        sprintf(buffer,"Sphere {}\n");
        UtilConcatVLA(&vla,&cc,buffer);
        UtilConcatVLA(&vla,&cc,"}\n\n");        
        break;
      case cPrimCylinder:
        break;
      case cPrimSausage:
        break;
      case cPrimTriangle:
        break;
      }
    }

    UtilConcatVLA(&vla,&cc,"}\n");
  }
  
  *vla_ptr=vla;
}
#endif


static int TriangleReverse(CPrimitive *p)
{
  float s1[3], s2[3], n0[3];

  subtract3f(p->v1,p->v2,s1);
  subtract3f(p->v3,p->v2,s2);
  cross_product3f(s1,s2,n0);
  
  if (dot_product3f(p->n0,n0) < 0)
    return 0;
  else
    return 1;
}

void RayRenderVRML2(CRay *I,int width,int height,
                    char **vla_ptr,float front,float back,
                    float fov, float angle,float z_corr)
{

  /* 

  From: pymol-users-admin@lists.sourceforge.net on behalf of Chris Want
  Sent: Tuesday, February 07, 2006 1:47 PM
  To: pymol-users@lists.sourceforge.net
  Subject: [PyMOL] VRML patch
  
  Hi Warren,
  
  I took your advice and modified the RayRenderVRML2() function to
  support triangles. I also threw out the sphere code that was already
  there and rewrote it (the code there was for VRML1, not
  VRML2). While I was at it, I also implemented export for cylinders
  and sausages.
  
  The code in the attached patch (diff-ed against cvs, and tested with
  two VRML2 readers) can be regarded as being in the public domain.
  
  Regards,
  Chris

  cwant_at_ualberta.ca

  */

  char *vla = *vla_ptr;
  int cc = 0; /* character count */
  OrthoLineType buffer;
  float mid[3]; /*, wid[3];*/
  float h_fov = cPI*(fov*width)/(180*height);
  RayExpandPrimitives(I);
  RayTransformFirst(I,0);
  RayComputeBox(I);
  /*
  mid[0] = (I->max_box[0] + I->min_box[0]) / 2.0;
  mid[1] = (I->max_box[1] + I->min_box[1]) / 2.0;
  mid[2] = (I->max_box[2] + I->min_box[2]) / 2.0;
  wid[0]  = (I->max_box[0] - I->min_box[0]);
  wid[1]  = (I->max_box[1] - I->min_box[1]);
  wid[2]  = (I->max_box[2] - I->min_box[2]);
  */

  copy3f(I->Pos,mid);
  UtilConcatVLA(&vla,&cc,
                "#VRML V2.0 utf8\n" /* WLD: most VRML2 readers req. utf8 */
                "\n");
  sprintf(buffer,
          "Viewpoint {\n"
          " position 0 0 %6.8f\n"
          " orientation 1 0 0 0\n"
          " description \"Z view\"\n"
          " fieldOfView %8.6f\n" /* WLD: use correct FOV */
          "}\n"
          /* WLD: only write the viewpoint which matches PyMOL
          "Viewpoint {\n"
          " position %6.8f 0 0\n"
          " orientation 0 1 0 1.570796\n"
          " description \"X view\"\n"
          "}\n"
          "Viewpoint {\n"
          " position 0 %6.8f 0\n"
          " orientation 0 -0.707106 -0.7071061 3.141592\n"
          " description \"Y view\"\n"
          "}\n"*/,
          -z_corr, /* *0.96646  for some reason, PyMOL and C4D cameras differ by about 3.5% ... */
          h_fov
          /*(wid[2] + wid[1]),
            (wid[0] + wid[1]),
            (wid[1] + wid[2])*/
          );
  UtilConcatVLA(&vla,&cc,buffer);
  {
    float light[3];
    float *lightv=SettingGetfv(I->G,cSetting_light);
    copy3f(lightv,light);
    normalize3f(light);
    sprintf(buffer,
            "DirectionalLight {\n"
            " direction %8.6f %8.6f %8.3f\n"
            "}\n",
            light[0],light[1],light[2]);
    UtilConcatVLA(&vla,&cc,buffer);
  }
  UtilConcatVLA(&vla,&cc,
                "NavigationInfo {\n"
                " headlight TRUE\n"
                " type \"EXAMINE\"\n"
                "}\n");
  { 
    int a, b;
    CPrimitive *prim;
    float *vert;
    int mesh_obj = false, mesh_start=0;

    CBasis *base = I->Basis+1;

    for(a=0;a<I->NPrimitive;a++) {
      prim = I->Primitive+a;
      vert = base->Vertex+3*(prim->vert);

      if(prim->type==cPrimTriangle) {
        if(!mesh_obj) {
          /* start mesh */
          mesh_start = a;
          UtilConcatVLA(&vla,&cc, 
                        "Shape {\n"
                        " appearance Appearance {\n"
                        "  material Material { diffuseColor 1.0 1.0 1.0 }\n"
                        " }\n"
                        " geometry IndexedFaceSet {\n"
                        "  coord Coordinate {\n"
                        "   point [\n");
          mesh_obj=true;
        }
      } else if(mesh_obj) {
        CPrimitive *cprim;
        int tri = 0;
        /* output connectivity */
        UtilConcatVLA(&vla,&cc, 
                      "   ]\n"
                      "  }\n"
                      "  coordIndex [\n");
        for(b=mesh_start;b<a;b++) {
          cprim = I->Primitive+b;
          if (TriangleReverse(cprim))
            sprintf(buffer,"%d %d %d -1,\n", tri, tri+2, tri+1);
          else
            sprintf(buffer,"%d %d %d -1,\n", tri, tri+1, tri+2);
          UtilConcatVLA(&vla,&cc,buffer);        
          tri+=3;
        }

        /* output vertex colors */
        UtilConcatVLA(&vla,&cc, 
                      "  ]\n"
                      "  colorPerVertex TRUE\n"
                      "  color Color {\n"
                      "   color [\n");
        for(b=mesh_start;b<a;b++) {
          cprim = I->Primitive+b;
          sprintf(buffer,
                  "%6.4f %6.4f %6.4f,\n"
                  "%6.4f %6.4f %6.4f,\n"
                  "%6.4f %6.4f %6.4f,\n", 
                  cprim->c1[0],cprim->c1[1],cprim->c1[2],
                  cprim->c2[0],cprim->c2[1],cprim->c2[2],
                  cprim->c3[0],cprim->c3[1],cprim->c3[2]);
          UtilConcatVLA(&vla,&cc,buffer);
        }

        /* output vertex normals */
        UtilConcatVLA(&vla,&cc, 
                      "  ] } \n"
                      "  normalPerVertex TRUE\n"
                      "  normal Normal {\n"
                      "   vector [\n");
        for(b=mesh_start;b<a;b++) {
          cprim = I->Primitive+b;
          {
            float *norm = base->Normal+3*base->Vert2Normal[cprim->vert];
            sprintf(buffer,
                    "%6.4f %6.4f %6.4f,\n"
                    "%6.4f %6.4f %6.4f,\n"
                  "%6.4f %6.4f %6.4f,\n", 
                    norm[3], norm[4], norm[5], /* transformed cprim->n1 */
                    norm[6], norm[7], norm[8], /* transformed cprim->n2 */
                    norm[9], norm[10], norm[11]); /* transformed cprim->n3 */
            UtilConcatVLA(&vla,&cc,buffer);
          }
        }
        UtilConcatVLA(&vla,&cc, 
                      "  ] }\n"
                      "  normalIndex [ \n");
        tri=0;
        for(b=mesh_start;b<a;b++) {
          cprim = I->Primitive+b;
          if (TriangleReverse(cprim))
            sprintf(buffer,"%d %d %d -1,\n", tri, tri+2, tri+1);
          else
            sprintf(buffer,"%d %d %d -1,\n", tri, tri+1, tri+2);
          UtilConcatVLA(&vla,&cc,buffer);
          tri+=3;
        }
        
        /* close mesh */
        UtilConcatVLA(&vla,&cc,
                      " ] \n"
                      " }\n"
                      "}\n");
        mesh_obj=false;
      }

      switch(prim->type) {
      case cPrimSphere:
        sprintf(buffer,
                "Transform {\n"
                " translation %8.6f %8.6f %8.6f\n"
                " children Shape {\n"
                "  geometry Sphere { radius %8.6f }\n"
                "  appearance Appearance {\n"
                "   material Material { diffuseColor %6.4f %6.4f %6.4f \n"
                "                       specularColor 0.8 0.8 0.8 \n"
                "                       shininess 8 }\n"
                "  }\n"
                " }\n"
                "}\n",        
                vert[0]-mid[0],
                vert[1]-mid[1],
                vert[2]-mid[2],
                prim->r1,
                prim->c1[0],prim->c1[1],prim->c1[2]);
        UtilConcatVLA(&vla,&cc,buffer);    
        break;
      case cPrimCylinder:
      case cPrimSausage:
        {
          float *d, vert2[3], axis[3], angle;
          OrthoLineType geometry;
          /* find the axis and angle that will rotate the y axis onto
           * the direction of the length of the cylinder
           */
          d=base->Normal+3*base->Vert2Normal[prim->vert];
          if ((d[0]*d[0] + d[2]*d[2]) < 0.000001) {
            /* parallel with y */
            axis[0] = 1.0;
            axis[1] = 0.0;
            axis[2] = 0.0;
            if (d[1] > 0) {
              angle = 0.0;
            }
            else {
              angle = cPI;
            }
          }
          else {
            axis[0] =  d[2];
            axis[1] =  0.0;
            axis[2] = -d[0];
            normalize3f(axis);
            angle = d[1];
            if (angle > 1.0) angle = 1.0;
            else if (angle < -1.0) angle = -1.0;
            angle = acos(angle);
          }
          /* vrml cylinders have origin in middle, not tip, that is why we
           * use prim->l1/2
           */
          scale3f(d,prim->l1/2,vert2);
          add3f(vert,vert2,vert2);
          if (prim->type==cPrimSausage) {
            OrthoLineType geom_add;
            sprintf(geometry,
                    "  Shape {\n"
                    "   geometry Cylinder {\n"
                    "    radius %8.6f\n"
                    "    height %8.6f\n"
                    "    bottom FALSE\n"
                    "    top    FALSE\n"
                    "   }\n"
                    "   appearance Appearance {\n"
                    "   material Material { diffuseColor %6.4f %6.4f %6.4f \n"
                    "                       specularColor 0.8 0.8 0.8 \n"
                    "                       shininess 8 }\n"
                    "   }\n"
                    "  }\n"
                    "  Transform {\n"
                    "   translation 0.0 %8.6f 0.0\n"
                    "   children Shape {\n"
                    "    geometry Sphere { radius %8.6f }\n"
                    "    appearance Appearance {\n"
                    "   material Material { diffuseColor %6.4f %6.4f %6.4f \n"
                    "                       specularColor 0.8 0.8 0.8 \n"
                    "                       shininess 8 }\n"
                    "    }\n"
                    "   }\n"
                    "  }\n",
                    prim->r1, prim->l1,
                    (prim->c1[0]+prim->c2[0])/2,
                    (prim->c1[1]+prim->c2[1])/2,
                    (prim->c1[2]+prim->c2[2])/2,
                    prim->l1/2, prim->r1,  
                    prim->c1[0],prim->c1[1],prim->c1[2]
                    );
            /* WLD: format string split to comply with ISO C89 standards */
            sprintf(geom_add,
                    "  Transform {\n"
                    "   translation 0.0 %8.6f 0.0\n"
                    "   children Shape {\n"
                    "    geometry Sphere { radius %8.6f }\n"
                    "    appearance Appearance {\n"
                    "   material Material { diffuseColor %6.4f %6.4f %6.4f \n"
                    "                       specularColor 0.8 0.8 0.8 \n"
                    "                       shininess 8 }\n"
                    "    }\n"
                    "   }\n"
                    "  }\n", 
                    -prim->l1/2, prim->r1,  
                    prim->c2[0],prim->c2[1],prim->c2[2]);
            strcat(geometry,geom_add);
          }
          else {
            sprintf(geometry,
                    "  Shape {\n"
                    "   geometry Cylinder {\n"
                    "    radius %8.6f\n"
                    "    height %8.6f\n"
                    "   }\n"
                    "   appearance Appearance {\n"
                    "   material Material { diffuseColor %6.4f %6.4f %6.4f \n"
                    "                       specularColor 0.8 0.8 0.8 \n"
                    "                       shininess 8 }\n"
                    "   }\n"
                    "  }\n",
                    prim->r1, prim->l1,
                    (prim->c1[0]+prim->c2[0])/2,
                    (prim->c1[1]+prim->c2[1])/2,
                    (prim->c1[2]+prim->c2[2])/2);
          }
          sprintf(buffer,
                  "Transform {\n"
                  " translation %8.6f %8.6f %8.6f\n"
                  " rotation %8.6f %8.6f %8.6f %8.6f\n"
                  " children [\n"
                  "%s"
                  " ]\n"
                  "}\n", 
                  vert2[0]-mid[0],
                  vert2[1]-mid[1],
                  vert2[2]-mid[2],
                  axis[0], axis[1], axis[2], angle,
                  geometry);
          UtilConcatVLA(&vla,&cc,buffer);
        }
        break;
      case cPrimTriangle:
        /* output coords. connectivity and vertex colors handled above/below */
        sprintf(buffer,
                "%8.6f %8.6f %8.6f,\n"
                "%8.6f %8.6f %8.6f,\n"
                "%8.6f %8.6f %8.6f,\n", 
                vert[0] - mid[0], vert[1] - mid[1], vert[2] - mid[2],
                vert[3] - mid[0], vert[4] - mid[1], vert[5] - mid[2],
                vert[6] - mid[0], vert[7] - mid[1], vert[8] - mid[2]);
        UtilConcatVLA(&vla,&cc,buffer);    
        break;
      }
    }

    if(mesh_obj) {
      CBasis *base = I->Basis+1;
      CPrimitive *cprim;
      int tri = 0;
      /* output connectivity */
      UtilConcatVLA(&vla,&cc, 
                    "   ]\n"
                    "  }\n"
                    "  coordIndex [\n");
      for(b=mesh_start;b<a;b++) {
        cprim = I->Primitive+b;
        if (TriangleReverse(cprim))
          sprintf(buffer,"%d %d %d -1,\n", tri, tri+2, tri+1);
        else
          sprintf(buffer,"%d %d %d -1,\n", tri, tri+1, tri+2);
        UtilConcatVLA(&vla,&cc,buffer);
        tri+=3;
      }

      /* output vertex colors */
      UtilConcatVLA(&vla,&cc,
                    "  ]\n"
                    "  colorPerVertex TRUE\n"
                    "  color Color {\n"
                    "   color [\n");
      for(b=mesh_start;b<a;b++) {
        cprim = I->Primitive+b;
        sprintf(buffer,
                "%6.4f %6.4f %6.4f,\n"
                "%6.4f %6.4f %6.4f,\n"
                "%6.4f %6.4f %6.4f,\n", 
                cprim->c1[0],cprim->c1[1],cprim->c1[2],
                cprim->c2[0],cprim->c2[1],cprim->c2[2],
                cprim->c3[0],cprim->c3[1],cprim->c3[2]);
        UtilConcatVLA(&vla,&cc,buffer);
      }

      /* output vertex normals */
      UtilConcatVLA(&vla,&cc, 
                    "  ] } \n"
                    "  normalPerVertex TRUE\n"
                    "  normal Normal {\n"
                    "   vector [\n");
      for(b=mesh_start;b<a;b++) {
        cprim = I->Primitive+b;
        {
          float *norm = base->Normal+3*base->Vert2Normal[cprim->vert];
          sprintf(buffer,
                  "%6.4f %6.4f %6.4f,\n"
                  "%6.4f %6.4f %6.4f,\n"
                  "%6.4f %6.4f %6.4f,\n", 
                  norm[3], norm[4], norm[5], /* transformed cprim->n1 */
                  norm[6], norm[7], norm[8], /* transformed cprim->n2 */
                  norm[9], norm[10], norm[11]); /* transformed cprim->n3 */
          UtilConcatVLA(&vla,&cc,buffer);
        }
      }
      UtilConcatVLA(&vla,&cc, 
                    "  ] }\n"
                    "  normalIndex [ \n");
      tri=0;
      for(b=mesh_start;b<a;b++) {
        cprim = I->Primitive+b;
        if (TriangleReverse(cprim))
          sprintf(buffer,"%d %d %d -1,\n", tri, tri+2, tri+1);
        else
          sprintf(buffer,"%d %d %d -1,\n", tri, tri+1, tri+2);
        UtilConcatVLA(&vla,&cc,buffer);
        tri+=3;
      }
    
      /* close mesh */
      UtilConcatVLA(&vla,&cc,
                    " ] \n"
                    " }\n"
                    "}\n");
      mesh_obj=false;
    }
  }
  
  *vla_ptr=vla;
}
/*========================================================================*/
void RayRenderObjMtl(CRay *I,int width,int height,char **objVLA_ptr,
                  char **mtlVLA_ptr,float front,float back,float fov,
                 float angle,float z_corr)
{
  char *objVLA = *objVLA_ptr; 
  char *mtlVLA = *mtlVLA_ptr; 
  int oc = 0; /* obj character count */
  /*  int mc = 0;*/ /* mtl character count */

  OrthoLineType buffer;
  
  RayExpandPrimitives(I);
  RayTransformFirst(I,0);

  { 
    int a;
    CPrimitive *prim;
    float *vert,*norm;
    int vc = 0;
    int nc = 0;
    CBasis *base = I->Basis+1;

    for(a=0;a<I->NPrimitive;a++) {
      prim = I->Primitive+a;
      vert = base->Vertex+3*(prim->vert);
      norm = base->Normal+3*base->Vert2Normal[prim->vert]+3;
      switch(prim->type) {
      case cPrimTriangle:
        sprintf(buffer,"v %8.6f %8.6f %8.6f\n",
                vert[0],vert[1],vert[2]-z_corr);
        UtilConcatVLA(&objVLA,&oc,buffer);
        sprintf(buffer,"v %8.6f %8.6f %8.6f\n",
                vert[3],vert[4],vert[5]-z_corr);
        UtilConcatVLA(&objVLA,&oc,buffer);
        sprintf(buffer,"v %8.6f %8.6f %8.6f\n",
                vert[6],vert[7],vert[8]-z_corr);
        UtilConcatVLA(&objVLA,&oc,buffer);
        sprintf(buffer,"vn %8.6f %8.6f %8.6f\n",
            norm[0],norm[1],norm[2]);
        UtilConcatVLA(&objVLA,&oc,buffer);
        sprintf(buffer,"vn %8.6f %8.6f %8.6f\n",
            norm[3],norm[4],norm[5]);
        UtilConcatVLA(&objVLA,&oc,buffer);
        sprintf(buffer,"vn %8.6f %8.6f %8.6f\n",
            norm[6],norm[7],norm[8]);
        UtilConcatVLA(&objVLA,&oc,buffer);
        if(TriangleReverse(prim)) {
          sprintf(buffer,"f %d//%d %d//%d %d//%d\n",
                  vc+1,nc+1,vc+3,nc+3,vc+2,nc+2);
        } else {
          sprintf(buffer,"f %d//%d %d//%d %d//%d\n",
                  vc+1,nc+1,vc+2,nc+2,vc+3,nc+3);
        }
        UtilConcatVLA(&objVLA,&oc,buffer);
        nc+=3;
        vc+=3;

        /*
        prim->c1[0],prim->c1[1],prim->c1[2])
        prim->c2[0],prim->c2[1],prim->c2[2],
        prim->c3[0],prim->c3[1],prim->c3[2]
        UtilConcatVLA(&vla,&oc,buffer);
        UtilConcatVLA(&vla,&oc,buffer);
        */

        break;
      case cPrimSphere:
      /*    sprintf(buffer,"v %8.6f %8.6f %8.6f\np %d\n",
                vert[0],vert[1],vert[2]-z_corr,vc+1);
            UtilConcatVLA(&objVLA,&oc,buffer);*/

      sprintf(buffer,"v %8.6f %8.6f %8.6f\n",
                vert[0],vert[1],vert[2]-z_corr);
        UtilConcatVLA(&objVLA,&oc,buffer);
      sprintf(buffer,"v %8.6f %8.6f %8.6f\n",
                vert[0],vert[1],vert[2]-z_corr);
        UtilConcatVLA(&objVLA,&oc,buffer);
      sprintf(buffer,"v %8.6f %8.6f %8.6f\n",
                vert[0],vert[1],vert[2]-z_corr);
        UtilConcatVLA(&objVLA,&oc,buffer);
        sprintf(buffer,"f %d %d %d\n",
                vc+1,vc+2,vc+3);
        UtilConcatVLA(&objVLA,&oc,buffer);
        vc+=3;
      break;
      }
    }
  }

  *objVLA_ptr = objVLA;
  *mtlVLA_ptr = mtlVLA;
}
/*========================================================================*/
void RayRenderPOV(CRay *I,int width,int height,char **headerVLA_ptr,
                  char **charVLA_ptr,float front,float back,float fov,
                  float angle,int antialias)
{
  int fogFlag=false;
  int fogRangeFlag=false;
  float fog;
  float *bkrd;
  float fog_start=0.0F;
  float gamma;
  float *d;
  CBasis *base;
  CPrimitive *prim;
  OrthoLineType buffer;
  float *vert,*norm;
  float vert2[3];
  int cc,hc;
  int a;
  int smooth_color_triangle;
  int mesh_obj = false;
  char *charVLA,*headerVLA;
  char transmit[64];
  float light[3],*lightv;
  float spec_power = SettingGet(I->G,cSetting_spec_power);
  if(spec_power<0.0F) {
    spec_power = SettingGet(I->G,cSetting_shininess);
  }
  spec_power/=4.0F;

  charVLA=*charVLA_ptr;
  headerVLA=*headerVLA_ptr;
  smooth_color_triangle=(int)SettingGet(I->G,cSetting_smooth_color_triangle);
  PRINTFB(I->G,FB_Ray,FB_Blather)
    " RayRenderPOV: w %d h %d f %8.3f b %8.3f\n",width,height,front,back
    ENDFB(I->G);
  if(Feedback(I->G,FB_Ray,FB_Blather)) {
    dump3f(I->Volume," RayRenderPOV: vol");
    dump3f(I->Volume+3," RayRenderPOV: vol");
  }
  cc=0;
  hc=0;
  gamma = SettingGet(I->G,cSetting_gamma);
  if(gamma>R_SMALL4)
    gamma=1.0F/gamma;
  else
    gamma=1.0F;

  lightv=SettingGetfv(I->G,cSetting_light);
  copy3f(lightv,light);

  fog = SettingGet(I->G,cSetting_ray_trace_fog);
  if(fog<0.0F)
    fog = SettingGet(I->G,cSetting_depth_cue);
  if(fog!=0.0F) {
    if(fog>1.0F) fog=1.0F;
    fogFlag=true;
    fog_start = SettingGet(I->G,cSetting_ray_trace_fog_start);
    if(fog_start<0.0F)
      fog_start = SettingGet(I->G,cSetting_fog_start);
    if(fog_start>1.0F)
      fog_start=1.0F;
    if(fog_start<0.0F)
      fog_start=0.0F;
    if(fog_start>R_SMALL4) {
      fogRangeFlag=true;
      if(fabs(fog_start-1.0F)<R_SMALL4) /* prevent div/0 */
        fogFlag=false;
    }
  }

  /* SETUP */
  
  if(antialias<0) 
    antialias = (int)SettingGet(I->G,cSetting_antialias);

  bkrd=SettingGetfv(I->G,cSetting_bg_rgb);

  RayExpandPrimitives(I);
  RayTransformFirst(I,0);

  PRINTFB(I->G,FB_Ray,FB_Details)
    " RayRenderPovRay: processed %i graphics primitives.\n",I->NPrimitive 
    ENDFB(I->G);
  base = I->Basis+1;

  if(!SettingGet(I->G,cSetting_ortho)) {
    sprintf(buffer,"camera {direction<0.0,0.0,%8.3f>\n location <0.0 , 0.0 , 0.0>\n right %12.10f*x up y \n }\n",
            -57.3F*cos(fov*cPI/(180*2.4))/fov,/* by trial and error */
            I->Range[0]/I->Range[1]);
  } else {
    sprintf(buffer,"camera {orthographic location <0.0 , 0.0 , %12.10f>\nlook_at  <0.0 , 0.0 , -1.0> right %12.10f*x up %12.10f*y}\n",
            front,-I->Range[0],I->Range[1]);
  }
  UtilConcatVLA(&headerVLA,&hc,buffer);

  {
    float ambient =           SettingGet(I->G,cSetting_ambient) + SettingGet(I->G,cSetting_direct);
    float reflect =           SettingGet(I->G,cSetting_reflect);

    if(ambient>0.5) ambient=0.5;

    reflect = 1.2F - 1.5F*ambient;

    sprintf(buffer,"#default { finish{phong %8.3f ambient %8.3f diffuse %8.3f phong_size %8.6f}}\n",
            SettingGet(I->G,cSetting_spec_reflect),
            ambient,
            reflect,
            spec_power);
    UtilConcatVLA(&headerVLA,&hc,buffer);
  }
  if(angle) {
    float temp[16];
    identity44f(temp);
    MatrixRotateC44f(temp,(float)-PI*angle/180,0.0F,1.0F,0.0F);
    MatrixTransformC44fAs33f3f(temp,light,light);
  }
  sprintf(buffer,"light_source{<%6.4f,%6.4f,%6.4f>  rgb<1.0,1.0,1.0>}\n",
          -light[0]*10000.0F,
          -light[1]*10000.0F,
          -light[2]*10000.0F-front
          );
  UtilConcatVLA(&headerVLA,&hc,buffer);

  {
    int opaque_back = SettingGetGlobal_i(I->G,cSetting_ray_opaque_background);
    if(opaque_back<0)
      opaque_back             = SettingGetGlobal_i(I->G,cSetting_opaque_background);      
    
    if(opaque_back) { /* drop a plane into the background for the background color */
      sprintf(buffer,
              "plane{z , %6.4f \n pigment{color rgb<%6.4f,%6.4f,%6.4f>}\n finish{phong 0 specular 0 diffuse 0 ambient 1.0}}\n"
              ,-back,bkrd[0],bkrd[1],bkrd[2]);
      UtilConcatVLA(&headerVLA,&hc,buffer);
    } 
  }
  
  for(a=0;a<I->NPrimitive;a++) {
    prim = I->Primitive+a;
    vert = base->Vertex+3*(prim->vert);
    if(prim->type==cPrimTriangle) {
      if(smooth_color_triangle)
        if(!mesh_obj) {
          sprintf(buffer,"mesh {\n");
          UtilConcatVLA(&charVLA,&cc,buffer);        
          mesh_obj=true;
        }
    } else if(mesh_obj) {
      sprintf(buffer," pigment{color rgb <1,1,1>}}");
      UtilConcatVLA(&charVLA,&cc,buffer);     
      mesh_obj=false;
    }
    switch(prim->type) {
       case cPrimSphere:
      sprintf(buffer,"sphere{<%12.10f,%12.10f,%12.10f>, %12.10f\n",
             vert[0],vert[1],vert[2],prim->r1);
      UtilConcatVLA(&charVLA,&cc,buffer);      
      sprintf(buffer,"pigment{color rgb<%6.4f,%6.4f,%6.4f>}}\n",
              prim->c1[0],prim->c1[1],prim->c1[2]);
      UtilConcatVLA(&charVLA,&cc,buffer);
            break;
       case cPrimCylinder:
      d=base->Normal+3*base->Vert2Normal[prim->vert];
      scale3f(d,prim->l1,vert2);
      add3f(vert,vert2,vert2);
      sprintf(buffer,"cylinder{<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>,\n %12.10f\n",
              vert[0],vert[1],vert[2],
              vert2[0],vert2[1],vert2[2],
              prim->r1);
      UtilConcatVLA(&charVLA,&cc,buffer);
      sprintf(buffer,"pigment{color rgb<%6.4f1,%6.4f,%6.4f>}}\n",
              (prim->c1[0]+prim->c2[0])/2,
              (prim->c1[1]+prim->c2[1])/2,
              (prim->c1[2]+prim->c2[2])/2);
      UtilConcatVLA(&charVLA,&cc,buffer);
            break;
    case cPrimSausage:
      d=base->Normal+3*base->Vert2Normal[prim->vert];
      scale3f(d,prim->l1,vert2);
      add3f(vert,vert2,vert2);
      sprintf(buffer,"cylinder{<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>,\n %12.10f\nopen\n",
              vert[0],vert[1],vert[2],
              vert2[0],vert2[1],vert2[2],
              prim->r1);
      UtilConcatVLA(&charVLA,&cc,buffer);
      sprintf(buffer,"pigment{color rgb<%6.4f1,%6.4f,%6.4f>}}\n",
              (prim->c1[0]+prim->c2[0])/2,
              (prim->c1[1]+prim->c2[1])/2,
              (prim->c1[2]+prim->c2[2])/2);
      UtilConcatVLA(&charVLA,&cc,buffer);

      sprintf(buffer,"sphere{<%12.10f,%12.10f,%12.10f>, %12.10f\n",
             vert[0],vert[1],vert[2],prim->r1);
      UtilConcatVLA(&charVLA,&cc,buffer);      
      sprintf(buffer,"pigment{color rgb<%6.4f1,%6.4f,%6.4f>}}\n",
              prim->c1[0],prim->c1[1],prim->c1[2]);
      UtilConcatVLA(&charVLA,&cc,buffer);

      sprintf(buffer,"sphere{<%12.10f,%12.10f,%12.10f>, %12.10f\n",
             vert2[0],vert2[1],vert2[2],prim->r1);
      UtilConcatVLA(&charVLA,&cc,buffer);      
      sprintf(buffer,"pigment{color rgb<%6.4f1,%6.4f,%6.4f>}}\n",
              prim->c2[0],prim->c2[1],prim->c2[2]);
      UtilConcatVLA(&charVLA,&cc,buffer);

      
            break;
       case cPrimTriangle:
      norm=base->Normal+3*base->Vert2Normal[prim->vert]+3;/* first normal is the average */      


      if(!TriangleDegenerate(vert,norm,vert+3,norm+3,vert+6,norm+6)) {

        if(smooth_color_triangle) {
          sprintf(buffer,"smooth_color_triangle{<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>,\n<%6.4f1,%6.4f,%6.4f>,\n<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>,\n<%6.4f1,%6.4f,%6.4f>,\n<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>,\n<%6.4f1,%6.4f,%6.4f> }\n",
                  vert[0],vert[1],vert[2],
                  norm[0],norm[1],norm[2],
                  prim->c1[0],prim->c1[1],prim->c1[2],
                  vert[3],vert[4],vert[5],
                  norm[3],norm[4],norm[5],
                  prim->c2[0],prim->c2[1],prim->c2[2],
                  vert[6],vert[7],vert[8],
                  norm[6],norm[7],norm[8],
                  prim->c3[0],prim->c3[1],prim->c3[2]
                  );
          UtilConcatVLA(&charVLA,&cc,buffer);
        } else {
#if 0
          
          sprintf(buffer,"smooth_triangle{<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>\n",
                  vert[0],vert[1],vert[2],
                  norm[0],norm[1],norm[2],
                  vert[3],vert[4],vert[5],
                  norm[3],norm[4],norm[5],
                  vert[6],vert[7],vert[8],
                  norm[6],norm[7],norm[8]
                  );
          UtilConcatVLA(&charVLA,&cc,buffer);
          if(prim->trans>R_SMALL4) 
            sprintf(transmit,"transmit %4.6f",prim->trans);
          else
            transmit[0]=0;
          if(equal3f(prim->c1,prim->c2)||equal3f(prim->c1,prim->c3)) {
            sprintf(buffer,"pigment{color rgb<%6.4f1,%6.4f,%6.4f> %s}}\n",
                    prim->c1[0],prim->c1[1],prim->c1[2],transmit);
          } else if(equal3f(prim->c2,prim->c3)) {
            sprintf(buffer,"pigment{color rgb<%6.4f1,%6.4f,%6.4f> %s}}\n",
                    prim->c2[0],prim->c2[1],prim->c2[2],transmit);
          } else {
            sprintf(buffer,"pigment{color rgb<%6.4f1,%6.4f,%6.4f> %s}}\n",
                    (prim->c1[0]+prim->c2[0]+prim->c3[0])/3,
                  (prim->c1[1]+prim->c2[1]+prim->c3[1])/3,
                    (prim->c1[2]+prim->c2[2]+prim->c3[2])/3,transmit);
          }
        UtilConcatVLA(&charVLA,&cc,buffer);
#else
        /* nowadays we use mesh2 to generate smooth_color_triangles */

        UtilConcatVLA(&charVLA,&cc,"mesh2 { ");
        sprintf(buffer,"vertex_vectors { 3, <%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>}\n normal_vectors { 3,\n<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>,\n<%12.10f,%12.10f,%12.10f>}\n",
                vert[0],vert[1],vert[2],
                vert[3],vert[4],vert[5],
                vert[6],vert[7],vert[8],
                norm[0],norm[1],norm[2],
                norm[3],norm[4],norm[5],
                norm[6],norm[7],norm[8]
                );
        UtilConcatVLA(&charVLA,&cc,buffer);

        if(prim->trans>R_SMALL4) 
          sprintf(transmit,"transmit %4.6f",prim->trans);
        else
          transmit[0]=0;

        sprintf(buffer,"texture_list { 3, ");
        UtilConcatVLA(&charVLA,&cc,buffer);

        sprintf(buffer, "texture { pigment{color rgb<%6.4f1,%6.4f,%6.4f> %s}}\n",
                prim->c1[0],prim->c1[1],prim->c1[2],transmit);
        UtilConcatVLA(&charVLA,&cc,buffer);

        sprintf(buffer, ",texture { pigment{color rgb<%6.4f1,%6.4f,%6.4f> %s}}\n",
                prim->c2[0],prim->c2[1],prim->c2[2],transmit);
        UtilConcatVLA(&charVLA,&cc,buffer);

        sprintf(buffer, ",texture { pigment{color rgb<%6.4f1,%6.4f,%6.4f> %s}} }\n",
                prim->c3[0],prim->c3[1],prim->c3[2],transmit);
        UtilConcatVLA(&charVLA,&cc,buffer);

        sprintf(buffer, "face_indices { 1, <0,1,2>, 0, 1, 2 } }\n");
        UtilConcatVLA(&charVLA,&cc,buffer);        
#endif

        }
      }
            break;
    }
  }
  
  if(mesh_obj) {
    sprintf(buffer," pigment{color rgb <1,1,1>}}");
    UtilConcatVLA(&charVLA,&cc,buffer);     
    mesh_obj=false;
  }
  *charVLA_ptr=charVLA;
  *headerVLA_ptr=headerVLA;
}

/*========================================================================*/
void RayProjectTriangle(CRay *I,RayInfo *r,float *light,float *v0,float *n0,float scale)
{
  register float w2;
  float d1[3],d2[3],d3[3];
  float p1[3],p2[3],p3[3];
  register int c=0;
  register const float _0 = 0.0F;
  register float *impact = r->impact;

  if(dot_product3f(light,n0-3)>=_0) c++;  
  else if(dot_product3f(light,n0)>=_0) c++;
  else if(dot_product3f(light,n0+3)>=_0) c++;
  else if(dot_product3f(light,n0+6)>=_0) c++;
  
  if(c) {

    w2 = 1.0F-(r->tri1+r->tri2);
    
    subtract3f(v0,impact,d1);
    subtract3f(v0+3,impact,d2);
    subtract3f(v0+6,impact,d3);
    project3f(d1,n0,p1);
    project3f(d2,n0+3,p2);
    project3f(d3,n0+6,p3);
    scale3f(p1,w2,d1);
    scale3f(p2,r->tri1,d2);
    scale3f(p3,r->tri2,d3);
    add3f(d1,d2,d2);
    add3f(d2,d3,d3);
    scale3f(d3,scale,d3);
    if(dot_product3f(r->surfnormal,d3)>=_0)
      add3f(d3,impact,impact);
  }
}
#ifndef _PYMOL_NOPY
static void RayHashSpawn(CRayHashThreadInfo *Thread,int n_thread,int n_total)
{
  int blocked;
  PyObject *info_list;
  int a,c,n=0;
  CRay *I = Thread->ray;
  PyMOLGlobals *G = I->G;

  blocked = PAutoBlock(G);

  PRINTFB(I->G,FB_Ray,FB_Blather)
    " Ray: filling voxels with %d threads...\n",n_thread
  ENDFB(I->G);
  while(n<n_total) {
    c = n;
    info_list = PyList_New(n_thread);
    for(a=0;a<n_thread;a++) {
      if((c+a)<n_total) {
        PyList_SetItem(info_list,a,PyCObject_FromVoidPtr(Thread+c+a,NULL));
      } else {
        PyList_SetItem(info_list,a,PConvAutoNone(NULL));
      }
      n++;
    }
    PXDecRef(PyObject_CallMethod(G->P_inst->cmd,"_ray_hash_spawn","OO",info_list,G->P_inst->cmd));
    Py_DECREF(info_list);
  }
  PAutoUnblock(G,blocked);
}
#endif

#ifndef _PYMOL_NOPY
static void RayAntiSpawn(CRayAntiThreadInfo *Thread,int n_thread)
{
  int blocked;
  PyObject *info_list;
  int a;
  CRay *I = Thread->ray;
  PyMOLGlobals *G = I->G;

  blocked = PAutoBlock(G);


  PRINTFB(I->G,FB_Ray,FB_Blather)
    " Ray: antialiasing with %d threads...\n",n_thread
  ENDFB(I->G);
  info_list = PyList_New(n_thread);
  for(a=0;a<n_thread;a++) {
    PyList_SetItem(info_list,a,PyCObject_FromVoidPtr(Thread+a,NULL));
  }
  PXDecRef(PyObject_CallMethod(G->P_inst->cmd,"_ray_anti_spawn","OO",info_list,G->P_inst->cmd));
  Py_DECREF(info_list);
  PAutoUnblock(G,blocked);
}
#endif

int RayHashThread(CRayHashThreadInfo *T)
{
  BasisMakeMap(T->basis,T->vert2prim,T->prim, T->n_prim, T->clipBox,T->phase,
               cCache_ray_map,T->perspective,T->front,T->size_hint);

  /* utilize a little extra wasted CPU time in thread 0 which computes the smaller map... */

  if(!T->phase) { 
    fill(T->image,T->background,T->bytes);
    RayComputeBox(T->ray);
  }
  return 1;
}
#ifndef _PYMOL_NOPY
static void RayTraceSpawn(CRayThreadInfo *Thread,int n_thread)
{
  int blocked;
  PyObject *info_list;
  int a;
  CRay *I=Thread->ray;
  PyMOLGlobals *G = I->G;

  blocked = PAutoBlock(G);

  PRINTFB(I->G,FB_Ray,FB_Blather)
    " Ray: rendering with %d threads...\n",n_thread
  ENDFB(I->G);
  info_list = PyList_New(n_thread);
  for(a=0;a<n_thread;a++) {
    PyList_SetItem(info_list,a,PyCObject_FromVoidPtr(Thread+a,NULL));
  }
  PXDecRef(PyObject_CallMethod(G->P_inst->cmd,"_ray_spawn","OO",info_list,G->P_inst->cmd));
  Py_DECREF(info_list);
  PAutoUnblock(G,blocked);
  
}
#endif

static int find_edge(unsigned int *ptr,float *depth, unsigned int width,
                     int threshold,int back)
{ /* can only be called for a pixel NOT on the edge */
  { /* color testing */
    register int compare0,compare1,compare2,compare3,compare4,compare5,compare6,compare7,compare8;
    {
      register int back_test, back_two = false;
      compare0 = (signed int)*(ptr);
      compare1 = (signed int)*(ptr-1);
      back_test = (compare0==back);
      compare2 = (signed int)*(ptr+1);
      back_two = back_two || ((compare1==back)==back_test);
      compare3 = (signed int)*(ptr-width);
      back_two = back_two || ((compare2==back)==back_test);
      compare4 = (signed int)*(ptr+width);
      back_two = back_two || ((compare3==back)==back_test);
      compare5 = (signed int)*(ptr-width-1);
      back_two = back_two || ((compare4==back)==back_test);
      compare6 = (signed int)*(ptr+width-1);
      back_two = back_two || ((compare5==back)==back_test);
      compare7 = (signed int)*(ptr-width+1);
      back_two = back_two || ((compare6==back)==back_test);
      compare8 = (signed int)*(ptr+width+1);
      back_two = back_two || ((compare7==back)==back_test);
      if(back_two) threshold = (threshold>>1); /* halve threshold for pixels that hit background */
    }
    
    {
      register int current;  
      register unsigned int shift = 0;
      register int sum1=0,sum2=3,sum3=0,sum4=0,sum5=0,sum6=0,sum7=0,sum8=0;
      int a;
      for(a=0;a<4;a++) {
        current = ((compare0>>shift)&0xFF);
        sum1 += abs(current - ((compare1>>shift)&0xFF));
        sum2 += abs(current - ((compare2>>shift)&0xFF));
        if(sum1>=threshold) return 1;
        sum3 += abs(current - ((compare3>>shift)&0xFF));
        if(sum2>=threshold) return 1;
        sum4 += abs(current - ((compare4>>shift)&0xFF));
        if(sum3>=threshold) return 1;
        sum5 += abs(current - ((compare5>>shift)&0xFF));
        if(sum4>=threshold) return 1;
        sum6 += abs(current - ((compare6>>shift)&0xFF));
        if(sum5>=threshold) return 1;
        sum7 += abs(current - ((compare7>>shift)&0xFF));
        if(sum6>=threshold) return 1;
        sum8 += abs(current - ((compare8>>shift)&0xFF));
        if(sum7>=threshold) return 1;
        if(sum8>=threshold) return 1;
        shift+=8;
      }
    }
  }
  if(depth) { /* depth testing */
    register float compare0,compare1,compare2,compare3,compare4,compare5,compare6,compare7,compare8;    
    register float dcutoff = threshold / 128.0F;
    
    compare1 = *(depth-1);
    compare0 = *(depth);
    compare2 = *(depth+1);
    if(fabs(compare0-compare1)>dcutoff) return 1;
    compare5 = *(depth-width-1);
    if(fabs(compare0-compare2)>dcutoff) return 1;
    compare3 = *(depth-width);
    if(fabs(compare0-compare5)>dcutoff) return 1;
    compare7 = *(depth-width+1);
    if(fabs(compare0-compare3)>dcutoff) return 1;
    compare6 = *(depth+width-1);
    if(fabs(compare0-compare7)>dcutoff) return 1;
    compare4 = *(depth+width);
    if(fabs(compare0-compare6)>dcutoff) return 1;
    compare8 = *(depth+width+1);
    if(fabs(compare0-compare4)>dcutoff) return 1;
    if(fabs(compare0-compare8)>dcutoff) return 1;
    /*    if(fabs(compare0-compare1)>0.001F)
      printf("%8.3f \n",compare0-compare1);
    if(fabs(compare0-compare2)>0.001F)
      printf("%8.3f \n",compare0-compare2);
    if(fabs(compare0-compare3)>0.001F)
      printf("%8.3f \n",compare0-compare3);
    if(fabs(compare0-compare4)>0.001F)
      printf("%8.3f \n",compare0-compare4);
    if(fabs(compare0-compare5)>0.001F)
      printf("%8.3f \n",compare0-compare5);
    if(fabs(compare0-compare6)>0.001F)
      printf("%8.3f \n",compare0-compare6);
    if(fabs(compare0-compare7)>0.001F)
      printf("%8.3f \n",compare0-compare7);
    if(fabs(compare0-compare8)>0.001F)
    printf("%8.3f \n",compare0-compare8);*/

  }
  return 0;
}

static void RayPrimGetColorRamped(PyMOLGlobals *G, float *matrix,RayInfo *r,float *fc)
{
  float fc1[3],fc2[3],fc3[3];
  register float *c1, *c2, *c3, w2;
  float back_pact[3];
  const float _0 = 0.0F, _1 = 1.0F, _01 = 0.1F;
  CPrimitive   *lprim   = r->prim;
  inverse_transformC44f3f(matrix,r->impact,back_pact);
  
  switch(lprim->type) {
  case cPrimTriangle:
    w2 = 1.0F - (r->tri1 + r->tri2);
    c1 = lprim->c1;
    if(c1[0]<=_0) {
      ColorGetRamped(G,(int)(c1[0]-_01),back_pact,fc1,-1);
      c1 = fc1;  
    }
    c2 = lprim->c2;
    if(c2[0]<=_0) {
      ColorGetRamped(G,(int)(c2[0]-_01),back_pact,fc2,-1);
      c2 = fc2;  
    }
    c3 = lprim->c3;
    if(c3[0]<=_0) {
      ColorGetRamped(G,(int)(c3[0]-_01),back_pact,fc3,-1);
      c3 = fc3;  
    }
    fc[0] = (c2[0]*r->tri1)+(c3[0]*r->tri2)+(c1[0]*w2);
    fc[1] = (c2[1]*r->tri1)+(c3[1]*r->tri2)+(c1[1]*w2);
    fc[2] = (c2[2]*r->tri1)+(c3[2]*r->tri2)+(c1[2]*w2);
    break;
  case cPrimSphere:
    c1 = lprim->c1;
    if(c1[0]<=_0) {
      ColorGetRamped(G,(int)(c1[0]-_01),back_pact,fc1,-1);
      c1 = fc1;  
    }
    copy3f(c1,fc);
    break;
  case cPrimEllipsoid:
    /* TO DO */
    break;
  case cPrimCylinder:
  case cPrimSausage:
    w2 = r->tri1;
    c1 = lprim->c1;
    if(c1[0]<=_0) {
      ColorGetRamped(G,(int)(c1[0]-_01),back_pact,fc1,-1);
      c1 = fc1;  
    }
    c2 = lprim->c2;
    if(c2[0]<=_0) {
      ColorGetRamped(G,(int)(c2[0]-_01),back_pact,fc2,-1);
      c2 = fc2;  
    }
    fc[0]=(c1[0]*(_1-w2))+(c2[0]*w2);
    fc[1]=(c1[1]*(_1-w2))+(c2[1]*w2);
    fc[2]=(c1[2]*(_1-w2))+(c2[2]*w2);
    break;
  default:
    fc[0] = _1;
    fc[1] = _1;
    fc[2] = _1;
    break;
  }
}

int RayTraceThread(CRayThreadInfo *T)
{
      CRay *I=T->ray;
      int x,y,yy;
      float excess=0.0F;
      float dotgle;
      float bright,direct_cmp,reflect_cmp,fc[4];
      float ambient,direct,lreflect,ft,ffact=0.0F,ffact1m;
      unsigned int cc0,cc1,cc2,cc3;
      int i;
      RayInfo r1,r2;
      int fogFlag=false;
      int fogRangeFlag=false;
      int opaque_back=0;
      int n_hit=0;
      int two_sided_lighting;
      float fog;
    float inter[3] = {0.0F,0.0F,0.0F};
      float fog_start=0.0F;
   /* float gamma,inp,sig=1.0F;*/
      float persist,persist_inv;
      float new_front;
      int pass;
      unsigned int last_pixel=0,*pixel;
      int exclude1,exclude2;
      float lit;
      int backface_cull;
      float project_triangle;
      float excl_trans;
      int shadows;
      int trans_shadows;
    int trans_mode;
      float first_excess;
      int pixel_flag;
      float ray_trans_spec, ray_lab_spec;
      float shadow_fudge;
    int label_shadow_mode;
      int interior_color;
      int interior_flag;
      int interior_shadows;
      int interior_wobble;
    int interior_mode;
   float interior_reflect;
      int wobble_save;
      float       settingPower, settingReflectPower,settingSpecPower,settingSpecReflect,settingSpecDirect;
    float       settingSpecDirectPower;
      float       invHgt, invFrontMinusBack, inv1minusFogStart,invWdth,invHgtRange;
      register float       invWdthRange,vol0;
      float       vol2;
      CBasis      *bp1,*bp2;
      int render_height;
      int offset=0;
   BasisCallRec BasisCall[MAX_BASIS];
   float border_offset;
   int edge_sampling = false;
   unsigned int edge_avg[4],edge_alpha_avg[4];
   int edge_cnt=0;
   float edge_base[2];
   float interior_normal[3];
   float edge_width = 0.35356F;
   float edge_height = 0.35356F;
   float trans_spec_cut,trans_spec_scale,trans_oblique,oblique_power;
   float direct_shade;
   float red_blend=0.0F;
   float blue_blend=0.0F;
   float green_blend=0.0F;
   float trans_cont;
   float pixel_base[3];
   float inv_trans_cont = 1.0F;
   float trans_cutoff,persist_cutoff;
   int trans_cont_flag = false;
   int blend_colors;
   int max_pass;
   float BasisFudge0,BasisFudge1;
   int perspective = T->perspective;
   float eye[3];
   float start[3],nudge[3],back_pact[3];
   float *depth = T->depth;
   float ray_scatter = SettingGetGlobal_f(I->G,cSetting_ray_scatter);
   const float shadow_decay = SettingGetGlobal_f(I->G,cSetting_ray_shadow_decay_factor);
   const float shadow_range = SettingGetGlobal_f(I->G,cSetting_ray_shadow_decay_range);
   const int clip_shadows = SettingGetGlobal_b(I->G,cSetting_ray_clip_shadows);
   const int spec_local = SettingGetGlobal_i(I->G,cSetting_ray_spec_local);
   float legacy = SettingGetGlobal_f(I->G,cSetting_ray_legacy_lighting);
   int spec_count = SettingGetGlobal_i(I->G,cSetting_spec_count);
   const float _0       = 0.0F;
   const float _1       = 1.0F;
   const float _p5            = 0.5F;
   const float _2       = 2.0F;
   const float _255     = 255.0F;
   const float _p499 = 0.499F;
   const float _persistLimit  = 0.0001F;
   float legacy_1m = _1 - legacy;
   int n_basis = I->NBasis;

   /*   MemoryDebugDump();
   printf("%d\n",sizeof(CPrimitive));
   */

   {
     float fudge = SettingGet(I->G,cSetting_ray_triangle_fudge);
     
     BasisFudge0 = 0.0F-fudge;
     BasisFudge1 = 1.0F+fudge;
   }
   if(spec_count<0) {
     spec_count = SettingGetGlobal_i(I->G,cSetting_light_count);
   }
   /* SETUP */
   
   /*  if(T->n_thread>1)
       printf(" Ray: Thread %d: Spawned.\n",T->phase+1);
   */
   
   interior_shadows     = SettingGetGlobal_i(I->G,cSetting_ray_interior_shadows);
   interior_wobble      = SettingGetGlobal_i(I->G,cSetting_ray_interior_texture);
   interior_color       = SettingGetGlobal_i(I->G,cSetting_ray_interior_color);
   interior_reflect  = 1.0F - SettingGet(I->G,cSetting_ray_interior_reflect);
   interior_mode = SettingGetGlobal_i(I->G,cSetting_ray_interior_mode);
   label_shadow_mode =  SettingGetGlobal_i(I->G,cSetting_label_shadow_mode);
   project_triangle     = SettingGet(I->G,cSetting_ray_improve_shadows);
   shadows                    = SettingGetGlobal_i(I->G,cSetting_ray_shadows);
   trans_shadows        = SettingGetGlobal_i(I->G,cSetting_ray_transparency_shadows);

   backface_cull        = SettingGetGlobal_i(I->G,cSetting_backface_cull);
   opaque_back                = SettingGetGlobal_i(I->G,cSetting_ray_opaque_background);
   if(opaque_back<0)
     opaque_back              = SettingGetGlobal_i(I->G,cSetting_opaque_background);      
   two_sided_lighting   = SettingGetGlobal_i(I->G,cSetting_two_sided_lighting);
   ray_trans_spec       = SettingGet(I->G,cSetting_ray_transparency_specular);
   ray_lab_spec         = SettingGet(I->G,cSetting_ray_label_specular);
   trans_cont        = SettingGetGlobal_f(I->G,cSetting_ray_transparency_contrast);
   trans_mode        = SettingGetGlobal_i(I->G,cSetting_transparency_mode);
   trans_oblique     = SettingGetGlobal_f(I->G,cSetting_ray_transparency_oblique);
   oblique_power     = SettingGetGlobal_f(I->G,cSetting_ray_transparency_oblique_power);
   trans_cutoff     = SettingGetGlobal_f(I->G,cSetting_ray_trace_trans_cutoff);
   persist_cutoff     = SettingGetGlobal_f(I->G,cSetting_ray_trace_persist_cutoff);

   if(trans_mode==1) two_sided_lighting = true;
   if(trans_cont>1.0F) {
     trans_cont_flag = true;
     inv_trans_cont = 1.0F/trans_cont;
   }
      ambient                       = T->ambient;
    /* divide up the reflected light component over all lights */
    {
      float reflect_scale = SceneGetReflectScaleValue(I->G,10);
      lreflect                = reflect_scale * (SettingGetGlobal_f(I->G,cSetting_reflect) - ray_scatter);
      if(lreflect<_0) lreflect=_0;
      ray_scatter = ray_scatter * reflect_scale;
    }
    direct                    = SettingGet(I->G,cSetting_direct);

      /* apply legacy adjustments */

      ambient*=(1.0F - legacy)+(legacy*(0.22F/0.12F));
      lreflect*=(1.0F - legacy)+(legacy*(0.72F/0.45F));
      direct*=(1.0F - legacy)+(legacy*(0.24F/0.45F));

      direct_shade      = SettingGet(I->G,cSetting_ray_direct_shade);
    trans_spec_cut = SettingGet(I->G,cSetting_ray_transparency_spec_cut);
    blend_colors    = SettingGetGlobal_i(I->G,cSetting_ray_blend_colors);
    max_pass = SettingGetGlobal_i(I->G,cSetting_ray_max_passes);
    if(blend_colors) {
      red_blend = SettingGet(I->G,cSetting_ray_blend_red);
     green_blend = SettingGet(I->G,cSetting_ray_blend_green);
     blue_blend = SettingGet(I->G,cSetting_ray_blend_blue);
   }

   if(trans_spec_cut<_1)
     trans_spec_scale = _1/(_1-trans_spec_cut);
   else
     trans_spec_scale = _0;

      /* COOP */
   settingPower         = SettingGet(I->G,cSetting_power);
   settingReflectPower  = SettingGet(I->G,cSetting_reflect_power);
   settingSpecPower     = SettingGet(I->G,cSetting_spec_power);
   if(settingSpecPower<0.0F) {
     settingSpecPower = SettingGet(I->G,cSetting_shininess);
   }
   
   {
     float spec_value = SettingGet(I->G,cSetting_specular);
     if(spec_value==1.0F) 
       spec_value = SettingGet(I->G,cSetting_specular_intensity);
     settingSpecReflect = SettingGet(I->G,cSetting_spec_reflect);
     if(settingSpecReflect<0.0F)
       settingSpecReflect = spec_value;     
     settingSpecReflect = SceneGetSpecularValue(I->G,settingSpecReflect,10);
     settingSpecDirect  = SettingGet(I->G,cSetting_spec_direct);
     if(settingSpecDirect<0.0F)
       settingSpecDirect = spec_value;
     settingSpecDirectPower   = SettingGet(I->G,cSetting_spec_direct_power);
     if(settingSpecDirectPower<0.0F)
       settingSpecDirectPower = settingSpecPower;
   }
   if(settingSpecReflect>1.0F) settingSpecReflect = 1.0F;
   if(SettingGet(I->G,cSetting_specular)<R_SMALL4) {
     settingSpecReflect = 0.0F;
   }
   
      if((interior_color!=-1)||(two_sided_lighting)||(trans_mode==1)||I->CheckInterior)
            backface_cull     = 0;

      shadow_fudge = SettingGet(I->G,cSetting_ray_shadow_fudge);
    
      inv1minusFogStart = _1;
      
      fog = SettingGet(I->G,cSetting_ray_trace_fog);
    if(fog<0.0F) {
      if(SettingGet(I->G,cSetting_depth_cue)) {
        fog = SettingGet(I->G,cSetting_fog);
      } else 
        fog = _0;
    }
   
      if(fog != _0) {
      if(fog>1.0F) fog=1.0F;
      fogFlag     = true;
      fog_start = SettingGet(I->G,cSetting_ray_trace_fog_start);
      if(fog_start<0.0F)
        fog_start = SettingGet(I->G,cSetting_fog_start);
      if(fog_start>1.0F)
        fog_start=1.0F;
      if(fog_start<0.0F)
        fog_start=0.0F;
      if(fog_start>R_SMALL4) {
        fogRangeFlag=true;
        if(fabs(fog_start-1.0F)<R_SMALL4) /* prevent div/0 */
          fogFlag=false;
      }
      inv1minusFogStart = _1 / (_1 - fog_start);
    }

    

      /* ray-trace */
      
   if(T->border) {
     invHgt                   = _1 / (float) (T->height-(3.0F+T->border));
     invWdth             = _1 / (float) (T->width-(3.0F+T->border));
   } else {

     invHgt                   = _1 / (float) (T->height);
     invWdth             = _1 / (float) (T->width);
   }

   if(perspective) {
     float height_range, width_range;

     zero3f(eye);

     /* subpixel-offsets for antialiasing naturally correspond to
        effective pixel sizes at the front of the visible slab...*/

     height_range = (T->front)*2*((float)tan((T->fov/2.0F)*PI/180.0F));
     width_range = height_range*(I->Range[0]/I->Range[1]);
     invWdthRange = invWdth * width_range;
     invHgtRange = invHgt * height_range;
     vol0 = eye[0] - width_range/2.0F;
     vol2 = eye[1] - height_range/2.0F;
   } else {
     invWdthRange        = invWdth * I->Range[0];
     invHgtRange         = invHgt * I->Range[1];
     vol0 = I->Volume[0];
     vol2 = I->Volume[2];
   }
      invFrontMinusBack = _1 / (T->front - T->back);

   edge_width *= invWdthRange;
   edge_height *= invHgtRange;

      bp1 = I->Basis + 1;
    if(I->NBasis>2) 
      bp2 = I->Basis + 2;
    else
      bp2 = NULL;

   render_height = T->y_stop - T->y_start;

   if(render_height) {
     offset = (T->phase * render_height/T->n_thread);
     offset = offset - (offset % T->n_thread) + T->phase;
   }
   if((interior_color!=-1)||I->CheckInterior) {


     if(interior_color!=-1)
       ColorGetEncoded(I->G,interior_color,inter);
     if(bp2) {
       interior_normal[0] = interior_reflect*bp2->LightNormal[0];
       interior_normal[1] = interior_reflect*bp2->LightNormal[1];
       interior_normal[2] = 1.0F+interior_reflect*bp2->LightNormal[2];
     } else {
       interior_normal[0] = 0.0;
       interior_normal[1] = 0.0;
       interior_normal[2] = 1.0F;
     }
     normalize3f(interior_normal);
   }
   
   r1.base[2]     = _0;

   BasisCall[0].Basis = I->Basis + 1;
   BasisCall[0].rr = &r1;
   BasisCall[0].vert2prim = I->Vert2Prim;
   BasisCall[0].prim = I->Primitive;
   BasisCall[0].shadow = false;
   BasisCall[0].back = T->back;
   BasisCall[0].trans_shadows = trans_shadows;
   BasisCall[0].nearest_shadow = (shadow_decay!=_0);
   BasisCall[0].check_interior = ((interior_color != -1) || I->CheckInterior);
   BasisCall[0].fudge0 = BasisFudge0;
   BasisCall[0].fudge1 = BasisFudge1;

   MapCacheInit(&BasisCall[0].cache,I->Basis[1].Map,T->phase,cCache_map_scene_cache);
   
   if(shadows&&(n_basis>2)) {
     int bc;
     for(bc=2;bc<n_basis;bc++) {
       BasisCall[bc].Basis = I->Basis + bc;
       BasisCall[bc].rr = &r2;
       BasisCall[bc].vert2prim = I->Vert2Prim;
       BasisCall[bc].prim = I->Primitive;
       BasisCall[bc].shadow = true;
       BasisCall[bc].front = _0;
       BasisCall[bc].back = _0;
       BasisCall[bc].excl_trans = _0;
       BasisCall[bc].trans_shadows = trans_shadows;
       BasisCall[bc].nearest_shadow =  (shadow_decay!=_0) || (clip_shadows);
       BasisCall[bc].check_interior = false;
       BasisCall[bc].fudge0 = BasisFudge0;
       BasisCall[bc].fudge1 = BasisFudge1;
       BasisCall[bc].label_shadow_mode = label_shadow_mode;
       MapCacheInit(&BasisCall[bc].cache,I->Basis[bc].Map,T->phase,cCache_map_shadow_cache);     
     }
   }
   
   if(T->border) {
     border_offset = -1.50F+T->border/2.0F;
   } else {
     border_offset = 0.0F;
   }
   for(yy = T->y_start; (yy < T->y_stop); yy++)
     {
       if(PyMOL_GetInterrupt(I->G->PyMOL,false))
         break;
       
       y = T->y_start + ((yy-T->y_start) + offset) % ( render_height); /* make sure threads write to different pages */
       
       if((!T->phase)&&!(yy & 0xF)) { /* don't slow down rendering too much */
         if(T->edging_cutoff) {
           if(T->edging) {
             OrthoBusyFast(I->G,(int)(2.5F*T->height/3 + 0.5F*y),4*T->height/3); 
           } else {
             OrthoBusyFast(I->G,(int)(T->height/3 + 0.5F*y),4*T->height/3); 
           }
         } else {
           OrthoBusyFast(I->G,T->height/3 + y,4*T->height/3); 
         }
       }
       pixel = T->image + (T->width * y) + T->x_start;
       
       if((y % T->n_thread) == T->phase) {      /* this is my scan line */
         pixel_base[1]  = ((y+0.5F+border_offset) * invHgtRange) + vol2;
           
         for(x = T->x_start; (x < T->x_stop); x++)
           {
               
             pixel_base[0]    = (((x+0.5F+border_offset)) * invWdthRange)  + vol0;
               
             while(1) {
               if(T->edging) {
                 if(!edge_sampling) {
                   if(x&&y&&(x<(T->width-1))&&(y<(T->height-1))) { /* not on the edge... */
                     if(find_edge(T->edging + (pixel - T->image),
                                  depth + (pixel - T->image),
                                  T->width, T->edging_cutoff,T->background)) {
                       register unsigned char *pixel_c = (unsigned char*)pixel;
                       register unsigned int c1,c2,c3,c4; 
                         
                       edge_cnt = 1;
                       edge_sampling = true;
                         
                       edge_avg[0] = (c1 = pixel_c[0]);
                       edge_avg[1] = (c2 = pixel_c[1]);
                       edge_avg[2] = (c3 = pixel_c[2]);
                       edge_avg[3] = (c4 = pixel_c[3]);
                         
                       edge_alpha_avg[0] = c1*c4;
                       edge_alpha_avg[1] = c2*c4;
                       edge_alpha_avg[2] = c3*c4;
                       edge_alpha_avg[3] = c4;
                         
                       edge_base[0]=pixel_base[0];
                       edge_base[1]=pixel_base[1];
                     }
                   }
                 }
                 if(edge_sampling) {
                   if(edge_cnt==5) {
                     /* done with edging, so store averaged value */

                     register unsigned char *pixel_c = (unsigned char*)pixel;
                     register unsigned int c1,c2,c3,c4; 

                     edge_sampling=false;
                     /* done with edging, so store averaged value */
                    
                     if(edge_alpha_avg[3]) {
                       c4 = edge_alpha_avg[3];
                       c1 = edge_alpha_avg[0] / c4;
                       c2 = edge_alpha_avg[1] / c4;
                       c3 = edge_alpha_avg[2] / c4;
                       c4 /= edge_cnt;
                     } else {
                       c1 = edge_avg[0]/edge_cnt;
                       c2 = edge_avg[1]/edge_cnt;
                       c3 = edge_avg[2]/edge_cnt;
                       c4 = edge_avg[3]/edge_cnt;
                     }
                     pixel_c[0] = c1;
                     pixel_c[1] = c2;
                     pixel_c[2] = c3;
                     pixel_c[3] = c4;

                     /* restore X,Y coordinates */
                     r1.base[0]=pixel_base[0];
                     r1.base[1]=pixel_base[1];

                   } else {
                     *pixel = T->background;
                     switch(edge_cnt) {
                     case 1:
                       r1.base[0] = edge_base[0]+edge_width;
                       r1.base[1] = edge_base[1]+edge_height;
                       break;
                     case 2:
                       r1.base[0] = edge_base[0]+edge_width;
                       r1.base[1] = edge_base[1]-edge_height;
                       break;
                     case 3:
                       r1.base[0] = edge_base[0]-edge_width;
                       r1.base[1] = edge_base[1]+edge_height;
                       break;
                     case 4:
                       r1.base[0] = edge_base[0]-edge_width;
                       r1.base[1] = edge_base[1]-edge_height;
                       break;
                     }

                   }
                 }
                 if(!edge_sampling) /* not oversampling this edge or already done... */
                   break;
               } else {
                 r1.base[0] = pixel_base[0];
                 r1.base[1] = pixel_base[1];
               }
              
               exclude1       = -1;
               exclude2     = -1;
               persist              = _1;
               first_excess   = _0;
               excl_trans           = _0;
               pass                 = 0;
               new_front            = T->front;

               if(perspective) {
                 r1.base[2] = -T->front;
                 r1.dir[0] = (r1.base[0] - eye[0]);
                 r1.dir[1] = (r1.base[1] - eye[1]);
                 r1.dir[2] = (r1.base[2] - eye[2]);
                 if(BasisCall[0].check_interior) {
                   start[0] = r1.base[0];
                   start[1] = r1.base[1];
                   start[2] = r1.base[2];
                 }
                 normalize3f(r1.dir);
                 {
                   register float scale = I->max_box[2]/r1.base[2];
                  
                   r1.skip[0] = r1.base[0]*scale;
                   r1.skip[1] = r1.base[1]*scale;
                   r1.skip[2] = I->max_box[2];
                 }

               }

               while((persist > _persistLimit) && (pass <= max_pass))
                 {
                   pixel_flag       = false;
                   BasisCall[0].except1 = exclude1;
                   BasisCall[0].except2 = exclude2;
                   BasisCall[0].front = new_front;
                   BasisCall[0].excl_trans = excl_trans;
                   BasisCall[0].interior_flag = false;

                   if(perspective) {
                     BasisCall[0].pass = pass;
                     if(pass) {
                       add3f(nudge,r1.base,r1.base);
                       copy3f(r1.base,r1.skip);
                     }
                     BasisCall[0].back_dist = -(T->back+r1.base[2])/r1.dir[2];
                     i = BasisHitPerspective( &BasisCall[0] );
                   } else {
                     i = BasisHitOrthoscopic( &BasisCall[0] );
                   }
                  
                   interior_flag = BasisCall[0].interior_flag;
                  
                   if(((i >= 0) || interior_flag) && (pass < max_pass))
                     {
                       pixel_flag         = true;
                       n_hit++;
                       if( ((r1.trans = r1.prim->trans) != _0 ) &&
                           trans_cont_flag ) {
                         r1.trans = (float)pow(r1.trans,inv_trans_cont);
                       }
                       if(interior_flag) {
                         copy3f(interior_normal,r1.surfnormal);
                         if(perspective) {
                           copy3f(start,r1.impact);
                           r1.dist = _0;
                         } else {
                           copy3f(r1.base,r1.impact);
                           r1.dist = T->front;
                           r1.impact[2]   -= T->front; 
                         }
                        
                         if(interior_wobble >= 0) {
                           wobble_save          = r1.prim->wobble; /* This is a no-no for multithreading! */
                           r1.prim->wobble      = interior_wobble;
                           
                           RayReflectAndTexture(I,&r1,perspective);
                           
                           r1.prim->wobble      = wobble_save;
                         } else
                           RayReflectAndTexture(I,&r1,perspective);
                        
                         dotgle = -r1.dotgle;
                         if((interior_color<0)&&(interior_color>cColorExtCutoff)) {
                           copy3f(r1.prim->ic,fc);
                         } else {
                           copy3f(inter,fc);
                         }
                       } else {
                         if(!perspective) 
                           new_front      = r1.dist;                        

                         switch(r1.prim->type) {
                         case cPrimTriangle:
                           
                           BasisGetTriangleNormal(bp1,&r1,i,fc,perspective);
                           r1.trans = (float)pow(r1.trans,inv_trans_cont);                         
 
                           if(r1.prim->ramped) {
                             RayPrimGetColorRamped(I->G, I->ModelView,&r1,fc);
                           }
                           if(bp2) {
                             RayProjectTriangle(I, &r1, bp2->LightNormal,
                                                bp1->Vertex+i*3,
                                                bp1->Normal+bp1->Vert2Normal[i]*3+3,
                                                project_triangle);
                           }
                          
                           RayReflectAndTexture(I,&r1,perspective);
                           if(perspective) {
                             BasisGetTriangleFlatDotglePerspective(bp1,&r1,i);
                           } else {
                             BasisGetTriangleFlatDotgle(bp1,&r1,i);
                           }
                           break;
                         case cPrimCharacter:
                           BasisGetTriangleNormal(bp1,&r1,i,fc,perspective);
                          
                           r1.trans = CharacterInterpolate(I->G,r1.prim->char_id,fc);
                          
                           RayReflectAndTexture(I,&r1,perspective);
                           BasisGetTriangleFlatDotgle(bp1,&r1,i);
                           break;

                         case cPrimEllipsoid:

                           BasisGetEllipsoidNormal(bp1,&r1,i,perspective);
                           RayReflectAndTexture(I,&r1,perspective);
                          
                           fc[0]=r1.prim->c1[0];
                           fc[1]=r1.prim->c1[1];
                           fc[2]=r1.prim->c1[2];
                           break;

                         default: /* sphere, cylinder, sausage, etc. */

                           /* must be a sphere (effectively speaking) */
                           
                           if(perspective) {
                             RayGetSphereNormalPerspective(I,&r1);
                           } else {
                             RayGetSphereNormal(I,&r1);
                           }

                           RayReflectAndTexture(I,&r1,perspective);
                          
                           if(r1.prim->ramped) {
                             RayPrimGetColorRamped(I->G, I->ModelView,&r1,fc);
                           } else if((r1.prim->type==cPrimCylinder) || (r1.prim->type==cPrimSausage)) {
                             ft = r1.tri1;
                             fc[0]=(r1.prim->c1[0]*(_1-ft))+(r1.prim->c2[0]*ft);
                             fc[1]=(r1.prim->c1[1]*(_1-ft))+(r1.prim->c2[1]*ft);
                             fc[2]=(r1.prim->c1[2]*(_1-ft))+(r1.prim->c2[2]*ft);
                           } else {
                             fc[0]=r1.prim->c1[0];
                             fc[1]=r1.prim->c1[1];
                             fc[2]=r1.prim->c1[2];
                           }
                           break;
                         }

                         if((trans_oblique!=_0)&&(r1.trans!=_0)) {
                           if(r1.surfnormal[2]>_0) {
                             float oblique_factor= r1.surfnormal[2];;
                             if(oblique_factor!=_1) {
                               if(r1.surfnormal[2]>_p5) {
                                 oblique_factor = (float)(_p5+_p5*(_1-pow((_1-oblique_factor)*_2,oblique_power)));
                               } else {
                                 oblique_factor = (float)(_p5*pow(oblique_factor*_2,oblique_power));
                               }
                             }
                             r1.trans *= (trans_oblique * oblique_factor) + (1.0F-trans_oblique);
                             if(r1.trans<0.06F) 
                               r1.trans=0.06F; /* don't allow transparent to become opaque */
                           }
                         }

                         dotgle=-r1.dotgle;
                        
                         if(r1.flat_dotgle < _0)
                           {
                             if((!two_sided_lighting) && (BasisCall[0].check_interior) && (interior_mode!=2))
                               {
                                 interior_flag        = true;
                                 copy3f(interior_normal,r1.surfnormal);
                                 if(perspective) {
                                   copy3f(start,r1.impact);                                    
                                   r1.dist = _0;
                                 } else {
                                   copy3f(r1.base,r1.impact);
                                   r1.impact[2]       -= T->front; 
                                   r1.dist                        = T->front;
                                 }
                                
                                 if(interior_wobble >= 0)
                                   {
                                     wobble_save            = r1.prim->wobble;
                                     r1.prim->wobble  = interior_wobble;
                                     RayReflectAndTexture(I,&r1,perspective);
                                     r1.prim->wobble  = wobble_save;
                                   }
                                 else
                                   RayReflectAndTexture(I,&r1,perspective);
                                
                                 dotgle   = -r1.dotgle;
                                 if((interior_color<0)&&(interior_color>cColorExtCutoff)) {
                                   copy3f(r1.prim->ic,fc);
                                 } else {
                                   copy3f(inter,fc);
                                 }

                               }
                           }
                        
                         if((dotgle < _0) && (!interior_flag))
                           {
                             if(two_sided_lighting) 
                               {
                                 dotgle   = -dotgle;
                                 invert3f(r1.surfnormal);
                               }
                             else 
                               dotgle     = _0;
                           }
                       }
                       
                       {
                         register double pow_dotgle;
                         register float pow_surfnormal2;
                         
                         if(settingPower!=_1) {
                           pow_dotgle = pow(dotgle, settingPower);
                           pow_surfnormal2 = (float)pow(r1.surfnormal[2], settingPower);
                         } else {
                           pow_dotgle = dotgle;
                           pow_surfnormal2 = r1.surfnormal[2];
                         }
                         direct_cmp = legacy_1m *  pow_surfnormal2 + /* new model */
                           legacy * ( (float) (dotgle + pow_dotgle) * _p5 ); /* legacy model */
                       }

                       reflect_cmp = _0;
                       if(settingSpecDirect!=_0) {
                         
#if 1
                         if(r1.surfnormal[2]>_0) {
                           excess   = (float)( pow(r1.surfnormal[2], settingSpecDirectPower) * settingSpecDirect);
                         } else {
                           excess =_0;
                         }
                                  
#else
                         float tmp[3];
                         tmp[0] = r1.dir[0];
                         tmp[1] = r1.dir[1];
                         tmp[2] = r2.dir[2]-_1;
                         dotgle     = -dot_product3f(r1.surfnormal,tmp);
                         if(dotgle < _0) dotgle=_0;                                                          
                         excess     = (float)( pow(dotgle, settingSpecDirectPower) * settingSpecDirect);
#endif

                       } else {
                         excess = _0;
                       }
                     
                       lit = _1;
                       if(n_basis<3) {
                         reflect_cmp = direct_cmp;
                       } else {
                         int bc;
                         CBasis *bp;
                         for(bc=2;bc<n_basis;bc++) {
                           lit = _1;
                           bp = I->Basis + bc;
                          
                           if(shadows && ((!interior_flag)||(interior_shadows)) &&
                              ((r1.prim->type != cPrimCharacter)||(label_shadow_mode&0x1))) {
                             matrix_transform33f3f(bp->Matrix,r1.impact,r2.base);
                             r2.base[2]-=shadow_fudge;
                             BasisCall[bc].except2 = -1;
                             BasisCall[bc].except1 = i; /* exclude current prim from shadow comp */
                             if(BasisHitShadow(&BasisCall[bc]) > -1) {
                               if( (!clip_shadows) || (bp->LightNormal[2]>=_0) || 
                                   ((T->front + r1.impact[2] - (r2.dist*bp->LightNormal[2]))<_0)) {
                                 lit = (float) pow(r2.trans, _p5);
                                 if((shadow_decay != _0) && (r2.dist>shadow_range)) {
                                   if(shadow_decay>0) {
                                     lit += ((_1-lit) * (_1 - _1 / exp((r2.dist-shadow_range) * shadow_decay)));
                                   } else {
                                     lit += ((_1-lit) * (_1 - _1 / pow(r2.dist/shadow_range,-shadow_decay)));
                                   }
                                 }
                               }
                             }
                           }
                          
                           if(lit>_0) {
                             {
                               register double pow_dotgle;
                               
                               dotgle     = -dot_product3f(r1.surfnormal,bp->LightNormal);
                               if(dotgle < _0) dotgle = _0;
                               
                               if(settingReflectPower!=_1)
                                 pow_dotgle = pow(dotgle, settingReflectPower);
                               else
                                 pow_dotgle = dotgle;
                               
                               reflect_cmp += legacy_1m * ((float)(lit * pow_dotgle)) + /* new model */
                                 legacy * ((float)(lit * (dotgle + pow_dotgle) * _p5 )); /* legacy model */
                             }

                             if(bc<(spec_count+2)) {

                               if(ray_scatter!=_0) /* scattered specular light */
                                 excess += ray_scatter * dotgle;

                               if(spec_local&&perspective) {
                                 /* slower, C4D-like local specular */
                                 float tmp[3];
                                 
                                 add3f(r1.surfnormal,r1.surfnormal,tmp);
                                 add3f(tmp,bp->LightNormal,tmp);
                                 normalize3f(tmp);
                                 dotgle   = -(dot_product3f(r1.dir,tmp))*1.004F;
                                 if(dotgle > _1) dotgle=_1;
                                 else if(dotgle < _0) dotgle=_0;
                                 dotgle = (float)( pow(dotgle, 0.29));
                               } else {
                                 dotgle   = -dot_product3f(r1.surfnormal,bp->SpecNormal); /* fast OpenGL-like global specular */
                               }
                               if(dotgle < _0) dotgle=_0;                                                          
                               if(r1.prim->type !=cPrimCharacter) {
                                 excess   += (float)( pow(dotgle, settingSpecPower) * settingSpecReflect * lit);
                               } else {
                                 excess   += (float)( pow(dotgle, settingSpecPower) * settingSpecReflect * lit * ray_lab_spec);
                               }
                             }
                           }
                         }
                       }
                       
                       if(fc[0]<=((float)cColorExtCutoff)) { /* ramped color */
                         inverse_transformC44f3f(I->ModelView,r1.impact,back_pact);
                         ColorGetRamped(I->G,(int)(fc[0]-0.1F),back_pact,fc,-1);
                       }
                      
                       bright = ambient +
                         (((_1-direct_shade)+direct_shade*lit) * direct*direct_cmp +
                          lreflect*reflect_cmp*(legacy_1m + legacy * direct_cmp)); /* blend legacy */
                       if(excess > _1) excess = _1;
                       if(bright > _1) bright = _1;
                       else if(bright < _0) bright = _0;
                      
                       /*                      bright *= (_1-excess);*/

                       fc[0] = (bright*fc[0]+excess);
                       fc[1] = (bright*fc[1]+excess);
                       fc[2] = (bright*fc[2]+excess);
                      
                       if(fogFlag) {
                         if(perspective) {
                           ffact = (T->front + r1.impact[2]) * invFrontMinusBack;
                         } else {
                           ffact = (T->front - r1.dist) * invFrontMinusBack;
                         }
                         if(fogRangeFlag)
                           ffact = (ffact - fog_start) * inv1minusFogStart;
                        
                         ffact*=fog;
                        
                         if(ffact<_0)     ffact = _0;
                         if(ffact>_1)     ffact = _1;
                        
                         ffact1m    = _1-ffact;
                        
                         if(opaque_back) {
                           fc[0]    = ffact*T->bkrd[0]+fc[0]*ffact1m;
                           fc[1]    = ffact*T->bkrd[1]+fc[1]*ffact1m;
                           fc[2]    = ffact*T->bkrd[2]+fc[2]*ffact1m;
                         } else {
                           fc[3] = ffact1m*(_1 - r1.trans);
                         }
                        
                         if(!pass) {
                           if(r1.trans<trans_spec_cut) {
                             first_excess = excess*ffact1m*ray_trans_spec;
                           } else {
                             first_excess = excess*ffact1m*ray_trans_spec*
                               trans_spec_scale*(_1 - r1.trans);
                           }
                         } else {
                           fc[0]+=first_excess; /* dubious? */
                           fc[1]+=first_excess;
                           fc[2]+=first_excess;
                         }
                       } else {
                         if(!pass) {
                           if(r1.trans<trans_spec_cut) {
                             first_excess = excess*ray_trans_spec;
                           } else {
                             first_excess = excess*ray_trans_spec*
                               trans_spec_scale*(_1 - r1.trans);
                           }
                         } else {
                           fc[0] += first_excess;
                           fc[1] += first_excess;
                           fc[2] += first_excess;
                         }
                         if(opaque_back) {
                           fc[3] = _1;
                         } else {
                           fc[3] = _1 - r1.trans;
                         }
                       }
                     }
                   else if(pass) {
                     /* hit nothing, and we're on on second or greater pass,
                        or we're on the last pass of a dead-end loop */
                     i=-1;
                     
                     fc[0] = first_excess+T->bkrd[0];
                     fc[1] = first_excess+T->bkrd[1];
                     fc[2] = first_excess+T->bkrd[2];
                     if(opaque_back) {
                       fc[3] = _1;
                     } else {
                       fc[3] = _0;
                     }
                     
                     ffact = 1.0F;
                     ffact1m = 0.0F;
                     
                     pixel_flag     = true;
                     if(trans_cont_flag)
                       persist = (float)pow(persist,trans_cont);
                     
                   }

                   if(pixel_flag)
                     {
                       /*
                         inp  = (fc[0]+fc[1]+fc[2]) * _inv3;
                         if(inp < R_SMALL4) 
                         sig = _1;
                         else
                         sig = (float)(pow(inp,gamma) / inp);
                      
                         cc0 = (uint)(sig * fc[0] * _255);
                         cc1 = (uint)(sig * fc[1] * _255);
                         cc2 = (uint)(sig * fc[2] * _255);
                       */

                       cc0 = (uint)(fc[0] * _255);
                       cc1 = (uint)(fc[1] * _255);
                       cc2 = (uint)(fc[2] * _255);

                       if(cc0 > 255) cc0 = 255;
                       if(cc1 > 255) cc1 = 255;
                       if(cc2 > 255) cc2 = 255;
                      
                       if(opaque_back) {
                         if(I->BigEndian) 
                           *pixel = T->fore_mask|(cc0<<24)|(cc1<<16)|(cc2<<8);
                         else
                           *pixel = T->fore_mask|(cc2<<16)|(cc1<<8)|cc0;
                       } else {
                         /* use alpha channel for fog with transparent backgrounds */
                         cc3  = (uint)(fc[3] * _255);
                         if(cc3 > 255) cc3 = 255;
                         
                         if(I->BigEndian)
                           *pixel = (cc0<<24)|(cc1<<16)|(cc2<<8)|cc3;
                         else
                           *pixel = (cc3<<24)|(cc2<<16)|(cc1<<8)|cc0;
                       }
                     }
                  
                   if(pass)   /* average all four channels */
                     {  
                       float mix_in;
                       if(i>=0) {
                         if(fogFlag) {
                           if(trans_cont_flag&&(ffact>_p5)) {
                             mix_in = 2*(persist*(_1-ffact)+((float)pow(persist,trans_cont)*(ffact-_p5)))
                               * (_1 - r1.trans*ffact);                            
                           } else {
                             mix_in = persist * (_1 - r1.trans*ffact);
                           }
                         } else {
                           mix_in = persist * (_1 - r1.trans);
                         }
                       } else {
                         mix_in = persist;
                       }

                       persist_inv = _1-mix_in;

                       if(!opaque_back) {
                         if(i<0) { /* hit nothing -- so don't blend */
                           fc[0] = (float)(0xFF&(last_pixel>>24));
                           fc[1] = (float)(0xFF&(last_pixel>>16));
                           fc[2] = (float)(0xFF&(last_pixel>>8));
                           fc[3] = (float)(0xFF&(last_pixel));
                           if(trans_cont_flag) { /* unless we are increasing contrast */
                             float m;
                             if(I->BigEndian) {
                               m = _1 - (float)(0xFF&(last_pixel))/_255;
                             } else {
                               m = _1 - (float)(0xFF&(last_pixel>>24))/_255;
                             }
                             m = _1 - (float)pow(m,trans_cont);
                             if(I->BigEndian) {
                               fc[3]      = m*_255 + _p499;
                             } else {
                               fc[0]      = m*_255 + _p499;
                             }
                           }
                         } else { /* hit something -- so keep blend and compute cumulative alpha*/
                          
                           fc[0]    = (0xFF&((*pixel)>>24)) * mix_in + (0xFF&(last_pixel>>24))*persist_inv;
                           fc[1]    = (0xFF&((*pixel)>>16)) * mix_in + (0xFF&(last_pixel>>16))*persist_inv;
                           fc[2]    = (0xFF&((*pixel)>>8))  * mix_in + (0xFF&(last_pixel>>8))*persist_inv;
                           fc[3]    = (0xFF&((*pixel)))     * mix_in + (0xFF&(last_pixel))*persist_inv;
                        
                           if(i>=0) { /* make sure opaque objects get opaque alpha*/
                             float o1,o2;
                             float m;
                            
                             if(I->BigEndian) {
                               o1 = (float)(0xFF&(last_pixel))/_255;
                               o2 = (float)(0xFF&(*pixel))/_255;
                             } else {
                               o1 = (float)(0xFF&(last_pixel>>24))/_255;
                               o2 = (float)(0xFF&((*pixel)>>24))/_255;
                             }
                            
                             if(o1<o2) { /* make sure o1 is largest opacity*/
                               m = o1;
                               o1 = o2;
                               o2 = m;
                             }
                             m = o1 + (1.0F - o1) * o2;
                             if(I->BigEndian) {
                               fc[3]      = m*_255 + _p499;
                             } else {
                               fc[0]      = m*_255 + _p499;
                             }
                           }
                         }
                       } else { /* opaque background, so just blend */
                         fc[0]      = (0xFF&((*pixel)>>24)) * mix_in + (0xFF&(last_pixel>>24))*persist_inv;
                         fc[1]      = (0xFF&((*pixel)>>16)) * mix_in + (0xFF&(last_pixel>>16))*persist_inv;
                         fc[2]      = (0xFF&((*pixel)>>8))  * mix_in + (0xFF&(last_pixel>>8))*persist_inv;
                         fc[3]      = (0xFF&((*pixel)))     * mix_in + (0xFF&(last_pixel))*persist_inv;
                       }
                      
                       cc0          = (uint)(fc[0]);
                       cc1          = (uint)(fc[1]);
                       cc2          = (uint)(fc[2]);
                       cc3          = (uint)(fc[3]);
                      
                       if(cc0 > 255) cc0  = 255;
                       if(cc1 > 255) cc1  = 255;
                       if(cc2 > 255) cc2  = 255;
                       if(cc3 > 255) cc3  = 255;
                      
                       *pixel = (cc0<<24)|(cc1<<16)|(cc2<<8)|cc3;
                      
                     }

                   if(depth&&(i>=0)&&
                      (r1.trans<trans_cutoff)&&
                      (persist>persist_cutoff)) {
                     depth[pixel - T->image] =(T->front + r1.impact[2]);
                   }
                  
                   if(i >= 0)
                     {
                       if(r1.prim->type == cPrimSausage) {  /* carry ray through the stick */
                         if(perspective) 
                           excl_trans = (2*r1.surfnormal[2]*r1.prim->r1/r1.dir[2]);                          
                         else
                           excl_trans = new_front+(2*r1.surfnormal[2]*r1.prim->r1);
                       }

                       if((!backface_cull)&&(trans_mode!=2))
                         persist    = persist * r1.trans;
                       else 
                         {
                           if((persist < 0.9999F) && (r1.trans>0.05F))  {
                             /* don't combine transparent surfaces */ 
                             *pixel = last_pixel;
                           } else {
                             persist      = persist * r1.trans;
                           }
                         }
                     }
                  

                   if( i < 0 ) {    /* nothing hit */
                     break;
                   } else {
                     if(perspective) {
                       if(r1.prim->type!=cPrimCharacter) {
                         float extend = r1.dist + 0.00001F;
                         scale3f(r1.dir, extend , nudge);
                       } else {
                         float extend = r1.dist;
                         scale3f(r1.dir, extend , nudge);
                       }
                     }
                     last_pixel     = *pixel;
                     exclude2 = exclude1;
                     exclude1       = i;
                     pass++;
                   }
                  
                 } /* end of ray while */

               if(blend_colors) {
                
                 float red_min = _0;
                 float green_min = _0;
                 float blue_min = _0;
                 float red_part;
                 float green_part;
                 float blue_part;

                 if(I->BigEndian) {
                   fc[0] = (float)(0xFF&(*pixel>>24));
                   fc[1] = (float)(0xFF&(*pixel>>16));
                   fc[2] = (float)(0xFF&(*pixel>>8));
                   cc3   =        (0xFF&(*pixel));
                 } else {
                   cc3   =        (0xFF&(*pixel>>24));
                   fc[2] = (float)(0xFF&(*pixel>>16));
                   fc[1] = (float)(0xFF&(*pixel>>8));
                   fc[0] = (float)(0xFF&(*pixel));
                 }

                 red_part = red_blend * fc[0];
                 green_part = green_blend * fc[1];
                 blue_part = blue_blend * fc[2];
                
                 red_min = (green_part>blue_part) ? green_part : blue_part;
                 green_min = (red_part>blue_part) ? red_part : blue_part;
                 blue_min = (green_part>red_part) ? green_part : red_part;
                
                 if(fc[0]<red_min) fc[0] = red_min;
                 if(fc[1]<green_min) fc[1] = green_min;
                 if(fc[2]<blue_min) fc[2] = blue_min;

                 cc0 = (uint)(fc[0]);
                 cc1 = (uint)(fc[1]);
                 cc2 = (uint)(fc[2]);
                
                 if(cc0 > 255) cc0 = 255;
                 if(cc1 > 255) cc1 = 255;
                 if(cc2 > 255) cc2 = 255;
                
                 if(I->BigEndian) 
                   *pixel = (cc0<<24)|(cc1<<16)|(cc2<<8)|cc3;
                 else
                   *pixel = (cc3<<24)|(cc2<<16)|(cc1<<8)|cc0;
               }
            
               if(!T->edging) break;
               /* if here, then we're edging...
                  so accumulate averages */
               { 
                
                 register unsigned char *pixel_c = (unsigned char*)pixel;
                 register unsigned int c1,c2,c3,c4; 
                
                 edge_avg[0] += (c1 = pixel_c[0]);
                 edge_avg[1] += (c2 = pixel_c[1]);
                 edge_avg[2] += (c3 = pixel_c[2]);
                 edge_avg[3] += (c4 = pixel_c[3]);
                
                 edge_alpha_avg[0] += c1*c4;
                 edge_alpha_avg[1] += c2*c4;
                 edge_alpha_avg[2] += c3*c4;
                 edge_alpha_avg[3] += c4;

                 edge_cnt++;
               }
              
             } /* end of edging while */
             pixel++;
           }      /* end of for */
         
       }    /* end of if */
            
     }      /* end of for */
      
   /*  if(T->n_thread>1) 
       printf(" Ray: Thread %d: Complete.\n",T->phase+1);*/
   MapCacheFree(&BasisCall[0].cache,T->phase,cCache_map_scene_cache);
      
   if(shadows&&(I->NBasis>2)) {
     int bc;
     for(bc=2;bc<I->NBasis;bc++) {
       MapCacheFree(&BasisCall[bc].cache,T->phase,cCache_map_shadow_cache);
     }
   }
      
   return (n_hit);
}

/* this is both an antialias and a slight blur */

/* for whatever reason, greatly GCC perfers a linear sequence of
   accumulates over a single large expression -- the difference is
   huge: over 10% !!! */


#define combine4by4(var,src,mask) { \
  var =  ((src)[0 ] & mask)   ; \
  var += ((src)[1 ] & mask)   ; \
  var += ((src)[2 ] & mask)   ; \
  var += ((src)[3 ] & mask)   ; \
  var += ((src)[4 ] & mask)   ; \
  var +=(((src)[5 ] & mask)*13) ; \
  var +=(((src)[6 ] & mask)*13) ; \
  var += ((src)[7 ] & mask)   ; \
  var += ((src)[8 ] & mask)    ; \
  var +=(((src)[9 ] & mask)*13) ; \
  var +=(((src)[10] & mask)*13) ; \
  var += ((src)[11] & mask)   ; \
  var += ((src)[12] & mask)   ; \
  var += ((src)[13] & mask)   ; \
  var += ((src)[14] & mask)   ; \
  var += ((src)[15] & mask)   ; \
  var = (var >> 6) & mask; \
}

#define combine5by5(var,src,mask) { \
  var =  ((src)[0 ] & mask)   ; \
  var += ((src)[1 ] & mask)   ; \
  var += ((src)[2 ] & mask)   ; \
  var += ((src)[3 ] & mask)   ; \
  var += ((src)[4 ] & mask)   ; \
  var += ((src)[5 ] & mask)   ; \
  var +=(((src)[6 ] & mask)*5); \
  var +=(((src)[7 ] & mask)*5); \
  var +=(((src)[8 ] & mask)*5); \
  var += ((src)[9 ] & mask)   ; \
  var += ((src)[10] & mask)   ; \
  var +=(((src)[11] & mask)*5); \
  var +=(((src)[12] & mask)*8); \
  var +=(((src)[13] & mask)*5); \
  var += ((src)[14] & mask)   ; \
  var += ((src)[15] & mask)   ; \
  var +=(((src)[16] & mask)*5); \
  var +=(((src)[17] & mask)*5); \
  var +=(((src)[18] & mask)*5); \
  var += ((src)[19] & mask)   ; \
  var += ((src)[20] & mask)   ; \
  var += ((src)[21] & mask)   ; \
  var += ((src)[22] & mask)   ; \
  var += ((src)[23] & mask)   ; \
  var += ((src)[24] & mask)   ; \
  var = (var >> 6) & mask; \
 }

#define combine6by6(var,src,mask) { \
  var =  ((src)[0 ] & mask)   ; \
  var += ((src)[1 ] & mask)   ; \
  var += ((src)[2 ] & mask)   ; \
  var += ((src)[3 ] & mask)   ; \
  var += ((src)[4 ] & mask)   ; \
  var += ((src)[5 ] & mask)   ; \
  var += ((src)[6 ] & mask)   ; \
  var +=(((src)[7 ] & mask)*5); \
  var +=(((src)[8 ] & mask)*7); \
  var +=(((src)[9 ] & mask)*7); \
  var +=(((src)[10] & mask)*5); \
  var += ((src)[11] & mask)   ; \
  var += ((src)[12] & mask)   ; \
  var +=(((src)[13] & mask)*7); \
  var +=(((src)[14] & mask)*8); \
  var +=(((src)[15] & mask)*8); \
  var +=(((src)[16] & mask)*7); \
  var += ((src)[17] & mask)   ; \
  var += ((src)[18] & mask)   ; \
  var +=(((src)[19] & mask)*7); \
  var +=(((src)[20] & mask)*8); \
  var +=(((src)[21] & mask)*8); \
  var +=(((src)[22] & mask)*7); \
  var += ((src)[23] & mask)   ; \
  var += ((src)[24] & mask)   ; \
  var +=(((src)[25] & mask)*5); \
  var +=(((src)[26] & mask)*7); \
  var +=(((src)[27] & mask)*7); \
  var +=(((src)[28] & mask)*5); \
  var += ((src)[29] & mask)   ; \
  var += ((src)[30] & mask)   ; \
  var += ((src)[31] & mask)   ; \
  var += ((src)[32] & mask)   ; \
  var += ((src)[33] & mask)   ; \
  var += ((src)[34] & mask)   ; \
  var += ((src)[35] & mask)   ; \
  var = (var >> 7) & mask; \
 }

#define m00FF 0x00FF
#define mFF00 0xFF00
#define mFFFF 0xFFFF

int RayAntiThread(CRayAntiThreadInfo *T)
{
      int         src_row_pixels;
      
      unsigned int *pSrc;
      unsigned int *pDst;
   /*   unsigned int m00FF=0x00FF,mFF00=0xFF00,mFFFF=0xFFFF;*/
      int width;
      int height;
      int x,y,yy;
      unsigned int *p;
      int offset = 0;
      CRay *I = T->ray;

      OrthoBusyFast(I->G,9,10);
      width = (T->width/T->mag) - 2;
      height = (T->height/T->mag) - 2;
      
      src_row_pixels    = T->width;

      offset = (T->phase * height)/T->n_thread;
      offset = offset - (offset % T->n_thread) + T->phase;

      for(yy = 0; yy< height; yy++ ) {
      y = (yy + offset) % height; /* make sure threads write to different pages */
      
      if((y % T->n_thread) == T->phase)   { /* this is my scan line */
        register unsigned long c1,c2,c3,c4,a;
        register unsigned char *c;
        
        pSrc      = T->image + src_row_pixels * (y*T->mag);
        pDst      = T->image_copy + width * y ; 
        switch(T->mag) {
        case 2: 
          {
            for(x = 0; x < width; x++) {
              
              c = (unsigned char*)( p = pSrc + (x * T->mag));
              c1 = c2 = c3 = c4 = a = 0;

              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              
              c = (unsigned char*)(p += src_row_pixels);
                
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*13); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*13); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              
              c = (unsigned char*)(p += src_row_pixels);

              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*13); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*13); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;

              c = (unsigned char*)(p += src_row_pixels);

              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;

              if(c4) {
                c1 /= c4;
                c2 /= c4;
                c3 /= c4;
              } else { /* compute straight RGB average */
                
                c = (unsigned char*)( p = pSrc + (x * T->mag));
                c1 = c2 = c3 = 0;

                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                
                c = (unsigned char*)(p += src_row_pixels);
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=13*c[0]; c2+=13*c[1]; c3+=13*c[2]; c+=4;
                c1+=13*c[0]; c2+=13*c[1]; c3+=13*c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                
                c = (unsigned char*)(p += src_row_pixels);
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=13*c[0]; c2+=13*c[1]; c3+=13*c[2]; c+=4;
                c1+=13*c[0]; c2+=13*c[1]; c3+=13*c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                
                c = (unsigned char*)(p += src_row_pixels);
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                
                c1 = c1>>6;
                c2 = c2>>6;
                c3 = c3>>6;
              }
              
              c = (unsigned char*)(pDst++);
              
              *(c++) = (unsigned char)c1;
              *(c++) = (unsigned char)c2;
              *(c++) = (unsigned char)c3;
              *(c++) = (unsigned char)(c4>>6);
            }
          }
          break;
        case 3:
          {
            for(x = 0; x < width; x++) {
              
              c = (unsigned char*)( p = pSrc + (x * T->mag));
              c1 = c2 = c3 = c4 = a = 0;
              
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              
              c = (unsigned char*)(p += src_row_pixels);
                
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              
              c = (unsigned char*)(p += src_row_pixels);

              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*8); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;

              c = (unsigned char*)(p += src_row_pixels);

              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              
              c = (unsigned char*)(p += src_row_pixels);

              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;

              if(c4) {
                c1 /= c4;
                c2 /= c4;
                c3 /= c4;
              } else { /* compute straight RGB average */
                
                c = (unsigned char*)( p = pSrc + (x * T->mag));
                c1 = c2 = c3 = 0;
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                
                c = (unsigned char*)(p += src_row_pixels);
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                
                c = (unsigned char*)(p += src_row_pixels);
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=8*c[0]; c2+=8*c[1]; c3+=8*c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                
                c = (unsigned char*)(p += src_row_pixels);
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;

                c = (unsigned char*)(p += src_row_pixels);
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                
                c1 = c1>>6;
                c2 = c2>>6;
                c3 = c3>>6;
              }
              
              c = (unsigned char*)(pDst++);
              
              *(c++) = (unsigned char)c1;
              *(c++) = (unsigned char)c2;
              *(c++) = (unsigned char)c3;
              *(c++) = (unsigned char)(c4>>6);
            }
          }
          break;
        case 4:
          {
            for(x = 0; x < width; x++) {
              
              c = (unsigned char*)( p = pSrc + (x * T->mag));
              c1 = c2 = c3 = c4 = a = 0;
              
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              
              c = (unsigned char*)(p += src_row_pixels);
                
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*7); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*7); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              
              c = (unsigned char*)(p += src_row_pixels);

              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*7); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*8); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*8); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*7); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              
              c = (unsigned char*)(p += src_row_pixels);

              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*7); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*8); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*8); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*7); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              
              c = (unsigned char*)(p += src_row_pixels);

              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*7); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*7); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]*5); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              
              c = (unsigned char*)(p += src_row_pixels);

              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;
              c4+=(a=c[3]); c1+=c[0]*a; c2+=c[1]*a; c3+=c[2]*a; c+=4;

              if(c4) {
                c1 /= c4;
                c2 /= c4;
                c3 /= c4;
              } else { /* compute straight RGB average */
                
                c = (unsigned char*)( p = pSrc + (x * T->mag));
                c1 = c2 = c3 = 0;
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                
                c = (unsigned char*)(p += src_row_pixels);
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=7*c[0]; c2+=7*c[1]; c3+=7*c[2]; c+=4;
                c1+=7*c[0]; c2+=7*c[1]; c3+=7*c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;


                c = (unsigned char*)(p += src_row_pixels);
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=7*c[0]; c2+=7*c[1]; c3+=7*c[2]; c+=4;
                c1+=8*c[0]; c2+=8*c[1]; c3+=8*c[2]; c+=4;
                c1+=8*c[0]; c2+=8*c[1]; c3+=8*c[2]; c+=4;
                c1+=7*c[0]; c2+=7*c[1]; c3+=7*c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;

                c = (unsigned char*)(p += src_row_pixels);
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=7*c[0]; c2+=7*c[1]; c3+=7*c[2]; c+=4;
                c1+=8*c[0]; c2+=8*c[1]; c3+=8*c[2]; c+=4;
                c1+=8*c[0]; c2+=8*c[1]; c3+=8*c[2]; c+=4;
                c1+=7*c[0]; c2+=7*c[1]; c3+=7*c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;


                c = (unsigned char*)(p += src_row_pixels);
                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=7*c[0]; c2+=7*c[1]; c3+=7*c[2]; c+=4;
                c1+=7*c[0]; c2+=7*c[1]; c3+=7*c[2]; c+=4;
                c1+=5*c[0]; c2+=5*c[1]; c3+=5*c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;

                c = (unsigned char*)(p += src_row_pixels);

                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;                
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                c1+=c[0]; c2+=c[1]; c3+=c[2]; c+=4;
                
                c1 = c1>>7;
                c2 = c2>>7;
                c3 = c3>>7;
              }
              
              c = (unsigned char*)(pDst++);
              
              *(c++) = (unsigned char)c1;
              *(c++) = (unsigned char)c2;
              *(c++) = (unsigned char)c3;
              *(c++) = (unsigned char)(c4>>7);
            }
          }
          break;

        }
      }
    }
    return 1;
}

#ifdef PROFILE_BASIS
extern int n_cells;
extern int n_prims;
extern int n_triangles;
extern int n_spheres;
extern int n_cylinders;
extern int n_sausages;
extern int n_skipped;
#endif

/*========================================================================*/
void RayRender(CRay *I,unsigned int *image,double timing,
               float angle,int antialias)
{
  int a;
  unsigned int *image_copy = NULL;
  unsigned int back_mask,fore_mask=0,trace_word=0;
  unsigned int background,buffer_size;
  int opaque_back=0;
  int n_hit=0;
  float *bkrd_ptr,bkrd[3];
  double now;
  int shadows;
  int n_thread;
  int mag=1;
  int oversample_cutoff;
  int perspective = SettingGetGlobal_i(I->G,cSetting_ray_orthoscopic);
  int n_light = SettingGetGlobal_i(I->G,cSetting_light_count);
  float ambient;
  float *depth = NULL;
  float front = I->Volume[4];
  float back = I->Volume[5];
  float fov  = I->Fov;
  float *pos = I->Pos;
  int width = I->Width;
  int height = I->Height;
  int trace_mode;
  const float _0 = 0.0F, _p499 = 0.499F;
  if(n_light>10) n_light = 10;
  
  if(perspective<0)
    perspective = SettingGetGlobal_b(I->G,cSetting_ortho);
  perspective = !perspective;

  VLACacheSize(I->G,I->Primitive,CPrimitive,I->NPrimitive,0,cCache_ray_primitive);
#ifdef PROFILE_BASIS
  n_cells = 0;
  n_prims = 0;
  n_triangles = 0;
  n_spheres = 0;
  n_cylinders = 0;
  n_sausages = 0;
  n_skipped = 0;
#endif

  n_thread  = SettingGetGlobal_i(I->G,cSetting_max_threads);
  if(n_thread<1)
    n_thread=1;
  if(n_thread>PYMOL_MAX_THREADS)
    n_thread = PYMOL_MAX_THREADS;
  opaque_back = SettingGetGlobal_i(I->G,cSetting_ray_opaque_background);
  if(opaque_back<0)
    opaque_back               = SettingGetGlobal_i(I->G,cSetting_opaque_background);      

  trace_mode  = SettingGetGlobal_i(I->G,cSetting_ray_trace_mode);

  shadows = SettingGetGlobal_i(I->G,cSetting_ray_shadows);

  oversample_cutoff = SettingGetGlobal_i(I->G,cSetting_ray_oversample_cutoff);

  if(antialias<0) {
    antialias = SettingGetGlobal_i(I->G,cSetting_antialias);
  }
  
  if(trace_mode && (antialias==1))
    antialias=2;
  else if(trace_mode && antialias)
    antialias++;

  if(antialias<0) antialias=0;
  if(antialias>4) antialias=4;

  if((!antialias) || trace_mode)
    oversample_cutoff = 0;

  mag = antialias;
  if(mag<1) mag=1;
  
  if(antialias>1) {
    width=(width+2)*mag;
    height=(height+2)*mag;
    image_copy = image;
    buffer_size = mag*mag*width*height;
    image = CacheAlloc(I->G,unsigned int,buffer_size,0,cCache_ray_antialias_buffer);
    ErrChkPtr(I->G,image);
  } else {
    buffer_size = width*height;
  }
  if(trace_mode) {
    depth = Calloc(float,width*height);
  } else if(oversample_cutoff) {
    depth = Calloc(float,width*height);
  }
  ambient = SettingGet(I->G,cSetting_ambient);
  
  bkrd_ptr=SettingGetfv(I->G,cSetting_bg_rgb);
  copy3f(bkrd_ptr,bkrd);
  { /* adjust bkrd and trace to offset the effect of gamma correction */
    float gamma = SettingGet(I->G,cSetting_gamma);
    {
      register float inp;
      register float sig;
      inp = (bkrd[0]+bkrd[1]+bkrd[2])/3.0F;
      if(inp < R_SMALL4) 
        sig = 1.0F;
      else
        sig = (float)(pow(inp,gamma))/inp;
      bkrd[0] *= sig;
      bkrd[1] *= sig;
      bkrd[2] *= sig;
      if(bkrd[0]>1.0F) bkrd[0] = 1.0F;
      if(bkrd[1]>1.0F) bkrd[1] = 1.0F;
      if(bkrd[2]>1.0F) bkrd[2] = 1.0F;
      
#if 0
      inp = ambient;
      if(inp < R_SMALL4) 
        sig = 1.0F;
      else
        sig = (float)(pow(inp,gamma))/inp;
      ambient *= sig;
      if(ambient>1.0f) ambient = 1.0F;
#endif
    }
    if(trace_mode) {
      register float inp;
      register float sig;
      int trace_color = SettingGetGlobal_color(I->G,cSetting_ray_trace_color);
      float trgb[3], *trgb_v = ColorGet(I->G,trace_color);
      copy3f(trgb_v,trgb);
      
      inp = (trgb[0]+trgb[1]+trgb[2])/3.0F;
      if(inp < R_SMALL4) 
        sig = 1.0F;
      else
        sig = (float)(pow(inp,gamma))/inp;
      trgb[0] *= sig;
      trgb[1] *= sig;
      trgb[2] *= sig;
      if(trgb[0]>1.0F) trgb[0] = 1.0F;
      if(trgb[1]>1.0F) trgb[1] = 1.0F;
      if(trgb[2]>1.0F) trgb[2] = 1.0F;
      
      if(I->BigEndian) {
        trace_word =
          ((0xFF& ((unsigned int)(trgb[0]*255+_p499))) <<24)|
          ((0xFF& ((unsigned int)(trgb[1]*255+_p499))) <<16)|
          ((0xFF& ((unsigned int)(trgb[2]*255+_p499))) <<8 )|
          0xFF;
      } else {
        trace_word = 
          0xFF000000 |
          ((0xFF& ((unsigned int)(trgb[2]*255+_p499))) <<16)|
          ((0xFF& ((unsigned int)(trgb[1]*255+_p499))) <<8)|
          ((0xFF& ((unsigned int)(trgb[0]*255+_p499))) );
      }
    }
  }
  if(opaque_back) {
    if(I->BigEndian)
      back_mask = 0x000000FF;
    else
      back_mask = 0xFF000000;
    fore_mask = back_mask;
  } else {
    if(I->BigEndian) {
      back_mask = 0x00000000;
    } else {
      back_mask = 0x00000000;
    }
  }
  if(I->BigEndian) {
     background = back_mask|
      ((0xFF& ((unsigned int)(bkrd[0]*255+_p499))) <<24)|
      ((0xFF& ((unsigned int)(bkrd[1]*255+_p499))) <<16)|
      ((0xFF& ((unsigned int)(bkrd[2]*255+_p499))) <<8 );
  } else {
    background = back_mask|
      ((0xFF& ((unsigned int)(bkrd[2]*255+_p499))) <<16)|
      ((0xFF& ((unsigned int)(bkrd[1]*255+_p499))) <<8)|
      ((0xFF& ((unsigned int)(bkrd[0]*255+_p499))) );
  }

  OrthoBusyFast(I->G,2,20);

  PRINTFB(I->G,FB_Ray,FB_Blather) 
    " RayNew: Background = %x %d %d %d\n",background,(int)(bkrd[0]*255),
    (int)(bkrd[1]*255),(int)(bkrd[2]*255)
    ENDFB(I->G);

  if(!I->NPrimitive) { /* nothing to render! */
    fill(image,background,width * (unsigned int)height);
  } else {
    
    if(I->PrimSizeCnt) {
      float factor = SettingGetGlobal_f(I->G,cSetting_ray_hint_camera);
      I->PrimSize = I->PrimSize/(I->PrimSizeCnt*factor);
      /*      printf("avg dist %8.7f\n",I->PrimSize);*/
    } else {
      I->PrimSize = 0.0F;
    }

    RayExpandPrimitives(I);
    RayTransformFirst(I,perspective);

    OrthoBusyFast(I->G,3,20);

    now = UtilGetSeconds(I->G)-timing;

       PRINTFB(I->G,FB_Ray,FB_Blather)
      " Ray: processed %i graphics primitives in %4.2f sec.\n",I->NPrimitive,now
      ENDFB(I->G);

     
    I->NBasis = n_light + 1; 
    if(I->NBasis>MAX_BASIS)
      I->NBasis = MAX_BASIS;
    if(I->NBasis<2) 
      I->NBasis = 2;
    { /* light sources */
      int bc;
      for(bc=2;bc<I->NBasis;bc++) {
        BasisInit(I->G,I->Basis+bc,bc);
        
        { /* setup light & rotate if necessary  */
          float light[3],*lightv;
          switch(bc) {
          default:
          case 2:
            lightv=SettingGetfv(I->G,cSetting_light);
            break;
          case 3:
            lightv=SettingGetfv(I->G,cSetting_light2);
            break;
          case 4:
            lightv=SettingGetfv(I->G,cSetting_light3);
            break;
          case 5:
            lightv=SettingGetfv(I->G,cSetting_light4);
            break;
          case 6:
            lightv=SettingGetfv(I->G,cSetting_light5);
            break;
          case 7:
            lightv=SettingGetfv(I->G,cSetting_light6);
            break;
          case 8:
            lightv=SettingGetfv(I->G,cSetting_light7);
            break;
          case 9:
            lightv=SettingGetfv(I->G,cSetting_light8);
            break;
          case 10:
            lightv=SettingGetfv(I->G,cSetting_light9);
            break;
          }
          copy3f(lightv,light);
          normalize3f(light);
          
          if(angle) {
            float temp[16];
            identity44f(temp);
            MatrixRotateC44f(temp,(float)-PI*angle/180,0.0F,1.0F,0.0F);
            MatrixTransformC44fAs33f3f(temp,light,light);
          }
          
          I->Basis[bc].LightNormal[0]=light[0];
          I->Basis[bc].LightNormal[1]=light[1];
          I->Basis[bc].LightNormal[2]=light[2];
          normalize3f(I->Basis[bc].LightNormal);
          
          {
            float spec_vector[3];
            copy3f(I->Basis[bc].LightNormal,spec_vector);
            spec_vector[2]--; 
            normalize3f(spec_vector);
            copy3f(spec_vector, I->Basis[bc].SpecNormal);
          }
        }
        
        if(shadows) { /* don't waste time on shadows unless needed */
          BasisSetupMatrix(I->Basis+bc);
          RayTransformBasis(I,I->Basis+bc,bc);
        }
      }
    }

    OrthoBusyFast(I->G,4,20);
#ifndef _PYMOL_NOPY
    if(shadows&&(n_thread>1)) { /* parallel execution */

      CRayHashThreadInfo *thread_info = Calloc(CRayHashThreadInfo,I->NBasis);
      
      /* rendering map */

      thread_info[0].basis = I->Basis+1;
      thread_info[0].vert2prim = I->Vert2Prim;
      thread_info[0].prim = I->Primitive;
      thread_info[0].n_prim = I->NPrimitive;
      thread_info[0].clipBox = I->Volume;
      thread_info[0].phase = 0;
      thread_info[0].perspective = perspective;
      thread_info[0].front = front;

      thread_info[0].image = image;
      thread_info[0].background = background;
      thread_info[0].bytes = width * (unsigned int)height;
      thread_info[0].ray = I; /* for compute box */
      thread_info[0].size_hint = I->PrimSize;
      /* shadow map */

      {
        int bc;
        float factor = SettingGetGlobal_f(I->G,cSetting_ray_hint_shadow);
        for(bc=2;bc<I->NBasis;bc++) {
          thread_info[bc-1].basis = I->Basis+bc;
          thread_info[bc-1].vert2prim = I->Vert2Prim;
          thread_info[bc-1].prim = I->Primitive;
          thread_info[bc-1].n_prim = I->NPrimitive;
          thread_info[bc-1].clipBox = NULL;
          thread_info[bc-1].phase = bc-1;
          thread_info[bc-1].perspective = false; 
          thread_info[bc-1].front = _0;
          /* allowing these maps to be more fine helps performance */
          thread_info[bc-1].size_hint = I->PrimSize*factor;
        }
      }

      /* NOTE that we're not limiting the number of threads in this phase
         under the assumption that it will usually just be a few threads */
      RayHashSpawn(thread_info, n_thread, I->NBasis-1); 
     
      FreeP(thread_info);
    } else
#endif
      { /* serial execution */
        BasisMakeMap(I->Basis+1,I->Vert2Prim,I->Primitive,I->NPrimitive,
                     I->Volume,0,cCache_ray_map,
                     perspective,front,I->PrimSize);
        if(shadows) {
          int bc;
          float factor = SettingGetGlobal_f(I->G,cSetting_ray_hint_shadow);
          for(bc=2;bc<I->NBasis;bc++) {
            BasisMakeMap(I->Basis+bc,I->Vert2Prim,I->Primitive,I->NPrimitive,
                         NULL,bc-1,cCache_ray_map,false,_0,I->PrimSize*factor);
          }
        }
        
        /* serial tasks which RayHashThread does in parallel mode using the first thread */
        
        fill(image,background,width * (unsigned int)height);
        RayComputeBox(I);
        
      }

    OrthoBusyFast(I->G,5,20);
    now = UtilGetSeconds(I->G)-timing;

#ifdef _MemoryDebug_ON
    if(Feedback(I->G,FB_Ray,FB_Debugging)) {
      MemoryDebugDump();
    }
    if(shadows) {
      PRINTFB(I->G,FB_Ray,FB_Blather)
        " Ray: voxels: [%4.2f:%dx%dx%d], [%4.2f:%dx%dx%d], %d MB, %4.2f sec.\n",
        I->Basis[1].Map->Div,   I->Basis[1].Map->Dim[0],
        I->Basis[1].Map->Dim[1],I->Basis[1].Map->Dim[2],
        I->Basis[2].Map->Div,   I->Basis[2].Map->Dim[0],
        I->Basis[2].Map->Dim[2],I->Basis[2].Map->Dim[2],
        (int)(MemoryDebugUsage()/(1024.0*1024)),
        now
        ENDFB(I->G);
    } else {
      PRINTFB(I->G,FB_Ray,FB_Blather)
        " Ray: voxels: [%4.2f:%dx%dx%d], %d MB, %4.2f sec.\n",
        I->Basis[1].Map->Div,   I->Basis[1].Map->Dim[0],
        I->Basis[1].Map->Dim[1],I->Basis[1].Map->Dim[2],
        (int)(MemoryDebugUsage()/(1024.0*1024)),
        now
        ENDFB(I->G);
    }
#else
    if(shadows) {
      PRINTFB(I->G,FB_Ray,FB_Blather)
        " Ray: voxels: [%4.2f:%dx%dx%d], [%4.2f:%dx%dx%d], %4.2f sec.\n",
        I->Basis[1].Map->Div,   I->Basis[1].Map->Dim[0],
        I->Basis[1].Map->Dim[1],I->Basis[1].Map->Dim[2],
        I->Basis[2].Map->Div,   I->Basis[2].Map->Dim[0],
        I->Basis[2].Map->Dim[2],I->Basis[2].Map->Dim[2],
        now
        ENDFB(I->G);
    } else {
      PRINTFB(I->G,FB_Ray,FB_Blather)
        " Ray: voxels: [%4.2f:%dx%dx%d], %4.2f sec.\n",
        I->Basis[1].Map->Div,   I->Basis[1].Map->Dim[0],
        I->Basis[1].Map->Dim[1],I->Basis[1].Map->Dim[2],
        now
        ENDFB(I->G);
    }

#endif
    /* IMAGING */
        
    {
            /* now spawn threads as needed */
      CRayThreadInfo *rt = Calloc(CRayThreadInfo,n_thread);
      
      int x_start=0,y_start=0;
      int x_stop=0,y_stop=0;
      float x_test=_0, y_test=_0;
      int x_pixel, y_pixel;

      if(perspective) {
        int c;

        if(I->min_box[2]>-front)
          I->min_box[2] = -front;
        if(I->max_box[2]>-front)
          I->max_box[2] = -front;

        for(c=0;c<4;c++) {
          switch(c) {
          case 0:
            x_test = -I->min_box[0]/I->min_box[2];
            y_test = -I->min_box[1]/I->min_box[2];
            break;
          case 1:
            x_test = -I->min_box[0]/I->max_box[2];
            y_test = -I->min_box[1]/I->max_box[2];
            break;
          case 2:
            x_test = -I->max_box[0]/I->min_box[2];
            y_test = -I->max_box[1]/I->min_box[2];
            break;
          case 3:
            x_test = -I->max_box[0]/I->max_box[2];
            y_test = -I->max_box[1]/I->max_box[2];
            break;
          }

          /* project onto back to get the effective range */

          x_pixel = (int)(width * (((x_test*I->Volume[5])-I->Volume[0])/I->Range[0])); 
          y_pixel = (int)(height * (((y_test*I->Volume[5])-I->Volume[2])/I->Range[1]));

          if(!c) {
            x_start = x_pixel;
            x_stop = x_pixel;
            y_start = y_pixel;
            y_stop = y_pixel;
          } else {
            if(x_start>x_pixel) x_start = x_pixel;
            if(x_stop<x_pixel) x_stop = x_pixel;
            if(y_start>y_pixel) y_start = y_pixel;
            if(y_stop<y_pixel) y_stop = y_pixel;
          }
        }
        x_start -=2;
        x_stop +=2;
        y_start -=2;
        y_stop +=2;

        /*
          x_start = 0; 
          y_start = 0;
          x_stop = width;
          y_stop = height;*/

      } else {
        x_start = (int)((width * (I->min_box[0] - I->Volume[0]))/I->Range[0]) - 2;
        x_stop  = (int)((width * (I->max_box[0] - I->Volume[0]))/I->Range[0]) + 2;
        
        y_stop = (int)((height * (I->max_box[1] - I->Volume[2]))/I->Range[1]) + 2;
        y_start  = (int)((height * (I->min_box[1] - I->Volume[2]))/I->Range[1]) - 2;
        
      }
      if(x_start<0) x_start = 0;
      if(y_start<0) y_start = 0;
      if(x_stop>width) x_stop = width;
      if(y_stop>height) y_stop = height;
      
      for(a=0;a<n_thread;a++) {
        rt[a].ray = I;
        rt[a].width = width;
        rt[a].height = height;
        rt[a].x_start = x_start;
        rt[a].x_stop = x_stop;
        rt[a].y_start = y_start;
        rt[a].y_stop = y_stop;
        rt[a].image = image;
        rt[a].border = mag-1;
        rt[a].front = front;
        rt[a].back = back;
        rt[a].fore_mask = fore_mask;
        rt[a].bkrd = bkrd;
        rt[a].ambient = ambient;
        rt[a].background = background;
        rt[a].phase = a;
        rt[a].n_thread = n_thread;
        rt[a].edging = NULL;
        rt[a].edging_cutoff = oversample_cutoff; /* info needed for busy indicator */
        rt[a].perspective = perspective;
        rt[a].fov = fov;
        rt[a].pos[2] = pos[2];
        rt[a].depth = depth;
      }
      
#ifndef _PYMOL_NOPY
            if(n_thread>1)
        RayTraceSpawn(rt,n_thread);
      else
#endif
        RayTraceThread(rt);

      if(oversample_cutoff) { /* perform edge oversampling, if requested */
        unsigned int *edging;

        edging = CacheAlloc(I->G,unsigned int,buffer_size,0,cCache_ray_edging_buffer);
        
        memcpy(edging,image,buffer_size*sizeof(unsigned int));

        for(a=0;a<n_thread;a++) {
          rt[a].edging = edging;
        }

#ifndef _PYMOL_NOPY
        if(n_thread>1)
          RayTraceSpawn(rt,n_thread);
        else 
#endif
          RayTraceThread(rt);

        CacheFreeP(I->G,edging,0,cCache_ray_edging_buffer,false);
      }
      FreeP(rt);
    }
  }
  
  if(depth && trace_mode) { 
    float *delta = Alloc(float,3*width*height);
    int x,y;
    {
      register int xc,yc;
      register float d,dzdx,dzdy, *p,*q,dd;
      p = depth;
      q = delta;
      for(y=0;y<height;y++) 
        for(x=0;x<width;x++) {
          dzdx = 0.0F;
          dzdy = 0.0F;
          xc = 0;
          yc = 0;
          d = *p;
          if(x) {
            dd = d - p[-1];
            dzdx += dd;
            xc++;
          }
          if(x<(width-1)) {
            dd = p[1] - d;
            if((!xc)||(fabs(dd)>fabs(dzdx)))
              dzdx = dd;
            xc = 1;
          }
          if(y) {
            dd = d - p[-width];
            dzdy+= dd;
            yc++;
          }
          if(y<(height-1)) {
            dd = p[width] - d;
            if((!yc)||(fabs(dd)>fabs(dzdy)))
              dzdy = dd;
            yc = 1;
          }
          p++;
          *(q++) = dzdx;
          *(q++) = dzdy;
          /*
          if(((x == y )||(x==width/2)||(y==height/2))&&((dzdx!=0.0F) || (dzdy!=0.0F))) {
            printf("%5d %5d : %8.3f %8.3f\n",y,x,dzdx,dzdy);
            }*/

          *(q++) = sqrt1f(dzdx*dzdx+dzdy*dzdy);
        }
    }
    
    {
      int i;
      {
        const float _1 = 1.0F;

        float invFrontMinusBack     = _1 / (front - back);
        float inv1minusFogStart     = _1;
        int fogFlag = false;
        float fog_start = 0.0F;
        int fogRangeFlag = false;
        float fog = SettingGet(I->G,cSetting_ray_trace_fog);
        if(fog<0.0F) {
          if(SettingGet(I->G,cSetting_depth_cue)) {
            fog = SettingGet(I->G,cSetting_fog);
          } else 
            fog = _0;
        }
        
        if(fog != _0) {
          if(fog>1.0F) fog=1.0F;
          fogFlag = true;
          fog_start = SettingGet(I->G,cSetting_ray_trace_fog_start);
          if(fog_start<0.0F)
            fog_start = SettingGet(I->G,cSetting_fog_start);
          if(fog_start>1.0F)
            fog_start=1.0F;
          if(fog_start<0.0F)
            fog_start=0.0F;
          if(fog_start>R_SMALL4) {
            fogRangeFlag=true;
            if(fabs(fog_start-1.0F)<R_SMALL4) /* prevent div/0 */
              fogFlag=false;
          }
          inv1minusFogStart   = _1 / (_1 - fog_start);
        }
        
        if(fogFlag) { /* make sure we have depth values at every potentially drawn pixel */
          float *tmp = Alloc(float,width*height);
          float dep;
          float *p,*q;
          int cnt;
          for(i=0;i<3;i++) { /* three passes required */
            p = depth;
            q = tmp;
            for(y=0;y<height;y++) 
              for(x=0;x<width;x++) {
                if(fabs(*p)<R_SMALL4) {
                  dep = 0.0F;
                  cnt = 0;
                  if(x) {
                    if(fabs(p[-1])>R_SMALL4) {
                      dep+=p[-1];
                      cnt++;
                    }
                  }
                  if(x<(width-1)) {
                    if(fabs(p[1])>R_SMALL4) {
                      dep+=p[1];
                      cnt++;
                    }
                  }
                  if(y) {
                    if(fabs(p[-width])>R_SMALL4) {
                      dep+=p[-width];
                      cnt++;
                    }
                  }
                  if(y<(height-1)) {
                    if(fabs(p[width])>R_SMALL4) {
                      dep+=p[width];
                      cnt++;
                    }
                  }
                  if(cnt) {
                    dep/=cnt;
                    *q = dep;
                  }
                } else {
                  *q = *p;
                }
                p++;
                q++;
              }
            p = tmp;
            tmp = depth;
            depth = p;
          }
          FreeP(tmp);
        }
        {
          unsigned int *q = image;
          float *p = delta;
          int dc = 0;
          int width3 = width*3;

          float slope_f = SettingGetGlobal_f(I->G,cSetting_ray_trace_slope_factor);
          float depth_f = SettingGetGlobal_f(I->G,cSetting_ray_trace_depth_factor);
          float disco_f = SettingGetGlobal_f(I->G,cSetting_ray_trace_disco_factor);
          register float diff,max_depth;
          register float dot,min_dot,max_slope, max_dz,max_pz;
          register float dx,dy,dz,px=0.0F,py=0.0F,pz=0.0F,ddx,ddy;
          const float _8 = 0.08F;
          const float _4 = 0.4F;
          const float _25 = 0.25F;
          const float _m25 = -0.25F;
          float disco_f_625 = disco_f*0.625F;
          float disco_f_5 = disco_f*0.5F;
          float disco_f_45 = disco_f*0.45F;

          {
            float gain = I->PixelRadius/(SettingGetGlobal_f(I->G,cSetting_ray_trace_gain)*I->Magnified);
            if(antialias)
              gain /= antialias;
            slope_f *= gain;
            depth_f *= gain;
            disco_f *= gain;
          }
          for(y=0;y<height;y++) 
            for(x=0;x<width;x++) {
              dc = 0;
              max_slope = 0.0F;   
              max_depth = 0.0F;
              min_dot = 1.0F;
              max_dz = 0.0F;
              max_pz = 0.0F;
              dx = p[0];
              dy = p[1];
              dz = p[2];
              for(i=0;i<8;i++) {
                switch(i) {
                case 0:
                  if(x) {
                    px = p[-3];
                    py = p[-2];
                    pz = p[-1];
                  }
                  break;
                case 1:
                  if(x<(width-1)) {
                    px = p[3];
                    py = p[4];
                    pz = p[5];
                  }
                  break;
                case 2:
                  if(y) { 
                    px = p[-width3];
                    py = p[-width3+1];
                    pz = p[-width3+2];
                  }
                  break;
                case 3:
                  if(y<(height-1)) {
                    px = p[width3];
                    py = p[width3+1];
                    pz = p[width3+2];
                  }
                  break;
                case 4:
                  if(x&&y) { 
                    px = p[-width3-3];
                    py = p[-width3-2];
                    pz = p[-width3-1];
                  }
                  break;
                case 5:
                  if(x&&(y<(height-1))) {
                    px = p[width3-3];
                    py = p[width3-2];
                    pz = p[width3-1];
                  }
                  break;
                case 6:
                  if(y&&(x<(width-1))) { 
                    px = p[-width3+3];
                    py = p[-width3+4];
                    pz = p[-width3+5];
                  }
                  break;
                case 7:
                  if((y<(height-1))&&(x<(width-1))) {
                    px = p[width3+3];
                    py = p[width3+4];
                    pz = p[width3+5];
                  }
                  break;
                }
                ddx = dx-px;
                ddy = dy-py;
                diff = ddx*ddx + ddy*ddy;
                if(max_depth<diff) max_depth = diff;
                if((dz>R_SMALL4)&&(pz>R_SMALL4)) { 
                  dot = (dx/dz)*(px/pz) + (dy/dz)*(py/pz);
                  if(dot<min_dot) {
                    min_dot = dot;
                    max_dz = dz;
                    max_pz = pz;
                  }
                }
                /*                if(dz>max_dz) max_dz = dz;
                                  if(pz>max_pz) max_pz = pz;*/
                diff = fabs(dz-pz);
                if(diff>max_slope)
                  max_slope = diff;
              }
              if((max_slope>(slope_f)) /* depth */
                 || (max_depth>(depth_f))/* slope */
                 /* gradient discontinuities -- could probably use more tuning...*/
                 || ((min_dot<_8)  &&((max_dz>disco_f)||(max_pz>disco_f)))
                 || ((min_dot<_4)  &&((max_dz>disco_f_625)||(max_pz>disco_f_625)))
                 || ((min_dot<_25) &&((max_dz>disco_f_5)||(max_pz>disco_f_5)))
                 || ((min_dot<_m25)&&((max_dz>disco_f_45)&&(max_pz>disco_f_45)))
                 ) {
                if(fogFlag) {
                  
                  float ffact =  depth[q-image] * invFrontMinusBack;
                  float ffact1m;
                  float fc[4];
                  unsigned int cc0,cc1,cc2,cc3;

                  if(fogRangeFlag)
                    ffact = (ffact - fog_start) * inv1minusFogStart;
                
                  ffact*=fog;
                
                  if(ffact<_0)      ffact = _0;
                  if(ffact>_1)      ffact = _1;
                
                  ffact1m     = _1-ffact;

                  if(opaque_back) {
                    fc[0]     = (0xFF&(background>>24)) * ffact + (0xFF&(trace_word>>24))*ffact1m;
                    fc[1]     = (0xFF&(background>>16)) * ffact + (0xFF&(trace_word>>16))*ffact1m;
                    fc[2]     = (0xFF&(background>>8))  * ffact + (0xFF&(trace_word>>8))*ffact1m;
                    fc[3]     = (0xFF&(background))    * ffact + (0xFF&(trace_word))*ffact1m;
                  } else { /* if non-opaque background, then use alpha to blend */
                    fc[1]     = (0xFF&(trace_word>>16));
                    fc[2]     = (0xFF&(trace_word>>8));
                    if(I->BigEndian) {
                      fc[0]   = (0xFF&(trace_word>>24));
                      fc[3]   = (0xFF&(background))    * ffact + (0xFF&(trace_word))*ffact1m;
                    } else {
                      fc[0]   = (0xFF&(background>>24)) * ffact + (0xFF&(trace_word>>24))*ffact1m;
                      fc[3]   = (0xFF&(trace_word));
                    }
                  }
                  cc0         = (uint)(fc[0]);
                  cc1         = (uint)(fc[1]);
                  cc2         = (uint)(fc[2]);
                  cc3         = (uint)(fc[3]);
                
                  if(cc0 > 255) cc0 = 255;
                  if(cc1 > 255) cc1 = 255;
                  if(cc2 > 255) cc2 = 255;
                  if(cc3 > 255) cc3 = 255;
                
                  *q = (cc0<<24)|(cc1<<16)|(cc2<<8)|cc3;
                } else {
                  *q = trace_word;
                }
              } else if(trace_mode==2) { /* only draw edge */
                *q = background;
              } else if(trace_mode==3) { /* quantize */
                *q = (*q & 0xC0C0C0C0);
                *q = *q | ((*q)>>2) | ((*q)>>4) | ((*q)>>6);
              }
              p+=3;
              q++;
            }
        }
      }
    }
    FreeP(delta);
  }
  
  if(antialias>1) {
    /* now spawn threads as needed */
    CRayAntiThreadInfo *rt = Calloc(CRayAntiThreadInfo, n_thread);
    
    for(a=0;a<n_thread;a++) {
      rt[a].width = width;
      rt[a].height = height;
      rt[a].image = image;
      rt[a].image_copy = image_copy;
      rt[a].phase = a;
      rt[a].mag = mag; /* fold magnification */
      rt[a].n_thread = n_thread;
      rt[a].ray = I;
    }
    
#ifndef _PYMOL_NOPY
    if(n_thread>1)
      RayAntiSpawn(rt,n_thread);
    else 
#endif
      RayAntiThread(rt);
    FreeP(rt);
    CacheFreeP(I->G,image,0,cCache_ray_antialias_buffer,false);
    image = image_copy;
  }

  PRINTFD(I->G,FB_Ray)
    " RayRender: n_hit %d\n",n_hit
    ENDFD;
#ifdef PROFILE_BASIS

  printf("int n_cells = %d;\nint n_prims = %d;\nint n_triangles = %8.3f;\nint n_spheres = %8.3f;\nint n_cylinders = %8.3f;\nint n_sausages = %8.3f;\nint n_skipped = %8.3f;\n",
         n_cells,
         n_prims,
         n_triangles/((float)n_cells),
         n_spheres/((float)n_cells),
         n_cylinders/((float)n_cells),
         n_sausages/((float)n_cells),
         n_skipped/((float)n_cells));
#endif
  FreeP(depth);
}

void RayRenderColorTable(CRay *I,int width,int height,int *image)
{
  int x,y;
  unsigned int r=0,g=0,b=0;
  unsigned int *pixel,mask,*p;

  if(I->BigEndian)
    mask = 0x000000FF;
  else
    mask = 0xFF000000;

  p=(unsigned int*)image; 
  for(x=0;x<width;x++)
    for(y=0;y<height;y++)
      *(p++)=mask;
  
  if((width>=512)&&(height>=512)) {
    
    
    for(y=0;y<512;y++) 
      for(x=0;x<512;x++)        
        {
          pixel = (unsigned int*) (image+((width)*y)+x);
          if(I->BigEndian) {
            *(pixel)=
              mask|(r<<24)|(g<<16)|(b<<8);
          } else {
            *(pixel)=
              mask|(b<<16)|(g<<8)|r;
          }
          b = b + 4;
          if(!(0xFF&b)) { 
            b=0;
            g=g+4;
            if(!(0xFF&g)) {           
              g=0;
              r=r+4;
            }
          }
        }
  }
}
/*========================================================================*/
void RayWobble(CRay *I,int mode,float *v)
{
  I->Wobble=mode;
  if(v) 
    copy3f(v,I->WobbleParam);
}
/*========================================================================*/
void RayTransparentf(CRay *I,float v)
{
  if(v>1.0F) v=1.0F;
  if(v<0.0F) v=0.0F;
  I->Trans=v;
}

void RayInteriorColor3fv(CRay *I,float *v,int passive)
{
  I->IntColor[0]=(*v++);
  I->IntColor[1]=(*v++);
  I->IntColor[2]=(*v++);
  if(!passive) 
    I->CheckInterior=true;
}
/*========================================================================*/
void RayColor3fv(CRay *I,float *v)
{
  I->CurColor[0]=(*v++);
  I->CurColor[1]=(*v++);
  I->CurColor[2]=(*v++);
}
/*========================================================================*/
void RaySphere3fv(CRay *I,float *v,float r)
{
  CPrimitive *p;

  float *vv;

  VLACacheCheck(I->G,I->Primitive,CPrimitive,I->NPrimitive,0,cCache_ray_primitive);
  p = I->Primitive+I->NPrimitive;
  
  p->type = cPrimSphere;
  p->r1=r;
  p->trans=I->Trans;
  p->wobble=I->Wobble;
  p->ramped=(I->CurColor[0]<0.0F);

  I->PrimSize += 2*r;
  I->PrimSizeCnt++;

  /* 
    copy3f(I->WobbleParam,p->wobble_param);*/
  vv=p->v1;
  (*vv++)=(*v++);
  (*vv++)=(*v++);
  (*vv++)=(*v++);

  vv=p->c1;
  v=I->CurColor;
  (*vv++)=(*v++);
  (*vv++)=(*v++);
  (*vv++)=(*v++);

  vv=p->ic;
  v=I->IntColor;
  (*vv++)=(*v++);
  (*vv++)=(*v++);
  (*vv++)=(*v++);

  if(I->TTTFlag) {
    transformTTT44f3f(I->TTT,p->v1,p->v1);
  }

  if(I->Context) {
    RayApplyContextToVertex(I,p->v1);
  }

  I->NPrimitive++;
}

void RayGetScaledAxes(CRay *I,float *xn,float *yn)
{
  float *v;
  float vt[3];
  float xn0[3] = {1.0F,0.0F,0.0F};
  float yn0[3] = {0.0F,1.0F,0.0F};
  float v_scale;

  v = TextGetPos(I->G);

  if(I->TTTFlag) {
    transformTTT44f3f(I->TTT,v, vt);
  } else {
    copy3f(v,vt);
  }

  v_scale =  RayGetScreenVertexScale(I,vt)/I->Sampling;

  RayApplyMatrixInverse33(1,(float3*)xn0,I->Rotation,(float3*)xn0);    
  RayApplyMatrixInverse33(1,(float3*)yn0,I->Rotation,(float3*)yn0);    
  
  scale3f(xn0,v_scale,xn); 
  scale3f(yn0,v_scale,yn); 
}

/*========================================================================*/
void RayCharacter(CRay *I,int char_id)
{
  CPrimitive *p;
  float *v;
  float vt[3];
  float *vv;
  float width,height;
  float v_scale;

  v = TextGetPos(I->G);
  VLACacheCheck(I->G,I->Primitive,CPrimitive,I->NPrimitive+1,0,cCache_ray_primitive);
  p = I->Primitive+I->NPrimitive;

  p->type = cPrimCharacter;
  p->trans=I->Trans;
  p->char_id = char_id;
  p->wobble=I->Wobble;
  p->ramped=0;
  /*
    copy3f(I->WobbleParam,p->wobble_param);*/

  vv=p->v1;
  (*vv++)=v[0];
  (*vv++)=v[1];
  (*vv++)=v[2];

  if(I->TTTFlag) {
    transformTTT44f3f(I->TTT,p->v1,p->v1);
  }
  /* what's the width of 1 screen window pixel at this point in space? */

  v_scale =  RayGetScreenVertexScale(I,p->v1)/I->Sampling;

  RayGetScreenVertexScale(I,vt);

  if(I->Context) {
    RayApplyContextToVertex(I,p->v1);
  }

  {
    float xn[3] = {1.0F,0.0F,0.0F};
    float yn[3] = {0.0F,1.0F,0.0F};
    float zn[3] = {0.0F,0.0F,1.0F};
    float sc[3];
    float scale;
    float xorig, yorig, advance;
    int width_i,height_i;
    CPrimitive *pp = p+1;

    RayApplyMatrixInverse33(1,(float3*)xn,I->Rotation,(float3*)xn);    
    RayApplyMatrixInverse33(1,(float3*)yn,I->Rotation,(float3*)yn);    
    RayApplyMatrixInverse33(1,(float3*)zn,I->Rotation,(float3*)zn);    

    CharacterGetGeometry(I->G,char_id,
                         &width_i, &height_i, 
                         &xorig, &yorig, &advance);
    width = (float)width_i;
    height = (float)height_i;

    scale = v_scale * advance;
    scale3f(xn,scale,vt); /* advance raster position in 3-space */
    add3f(v,vt,vt);
    TextSetPos(I->G,vt);

    /* position the pixmap relative to raster position */

  
    /*    scale = ((-xorig)-0.5F)*I->PixelRadius;*/
    scale = ((-xorig)-0.0F)*v_scale;
    scale3f(xn,scale,sc);
    add3f(sc,p->v1,p->v1);
         
    scale = ((-yorig)-0.0F)*v_scale;
    scale3f(yn,scale,sc);
    add3f(sc,p->v1,p->v1);
    

    scale = v_scale*width;
    scale3f(xn,scale,xn);
    scale = v_scale*height;
    scale3f(yn,scale,yn);

    copy3f(zn,p->n0);
    copy3f(zn,p->n1);
    copy3f(zn,p->n2);
    copy3f(zn,p->n3);

    *(pp)=(*p);

    /* define coordinates of first triangle */

    add3f(p->v1,xn,p->v2);
    add3f(p->v1,yn,p->v3);

    I->PrimSize += 2*(diff3f(p->v1,p->v2) + diff3f(p->v1,p->v3) + diff3f(p->v2,p->v3));
    I->PrimSizeCnt += 6;
    
    /* encode characters coordinates in the colors  */

    zero3f(p->c1);
    set3f(p->c2,width,0.0F,0.0F);
    set3f(p->c3,0.0F,height,0.0F);

    /* define coordinates of second triangle */

    add3f(yn,xn,pp->v1);
    add3f(p->v1,pp->v1,pp->v1);
    add3f(p->v1,yn,pp->v2);
    add3f(p->v1,xn,pp->v3);

    {
      float *v,*vv;
      vv=p->ic;
      v=I->IntColor;
      (*vv++)=(*v++);
      (*vv++)=(*v++);
      (*vv++)=(*v++);
      vv=pp->ic;
      v=I->IntColor;
      (*vv++)=(*v++);
      (*vv++)=(*v++);
      (*vv++)=(*v++);
    }

    /* encode integral character coordinates into the vertex colors  */

    set3f(pp->c1,width,height,0.0F);
    set3f(pp->c2,0.0F,height,0.0F);
    set3f(pp->c3,width,0.0F,0.0F);

  }

  I->NPrimitive+=2;
}
/*========================================================================*/
void RayCylinder3fv(CRay *I,float *v1,float *v2,float r,float *c1,float *c2)
{
  CPrimitive *p;

  float *vv;

  VLACacheCheck(I->G,I->Primitive,CPrimitive,I->NPrimitive,0,cCache_ray_primitive);
  p = I->Primitive+I->NPrimitive;

  p->type = cPrimCylinder;
  p->r1=r;
  p->trans=I->Trans;
  p->cap1=cCylCapFlat;
  p->cap2=cCylCapFlat;
  p->wobble=I->Wobble;
  p->ramped = ((c1[0]<0.0F)||(c2[0]<0.0F));
  /* 
 copy3f(I->WobbleParam,p->wobble_param);*/

  vv=p->v1;
  (*vv++)=(*v1++);
  (*vv++)=(*v1++);
  (*vv++)=(*v1++);
  vv=p->v2;
  (*vv++)=(*v2++);
  (*vv++)=(*v2++);
  (*vv++)=(*v2++);

  I->PrimSize += diff3f(p->v1,p->v2) + 2*r;
  I->PrimSizeCnt++;

  if(I->TTTFlag) {
    transformTTT44f3f(I->TTT,p->v1,p->v1);
    transformTTT44f3f(I->TTT,p->v2,p->v2);
  }

  if(I->Context) {
    RayApplyContextToVertex(I,p->v1);
    RayApplyContextToVertex(I,p->v2);
  }

  vv=p->c1;
  (*vv++)=(*c1++);
  (*vv++)=(*c1++);
  (*vv++)=(*c1++);
  vv=p->c2;
  (*vv++)=(*c2++);
  (*vv++)=(*c2++);
  (*vv++)=(*c2++);

  {
    float *v;
    vv=p->ic;
    v=I->IntColor;
    (*vv++)=(*v++);
    (*vv++)=(*v++);
    (*vv++)=(*v++);
  }

  I->NPrimitive++;
}
/*========================================================================*/
void RayCustomCylinder3fv(CRay *I,float *v1,float *v2,float r,
                          float *c1,float *c2,int cap1,int cap2)
{
  CPrimitive *p;

  float *vv;

  VLACacheCheck(I->G,I->Primitive,CPrimitive,I->NPrimitive,0,cCache_ray_primitive);
  p = I->Primitive+I->NPrimitive;

  p->type = cPrimCylinder;
  p->r1=r;
  p->trans=I->Trans;
  p->cap1=cap1;
  p->cap2=cap2;
  p->wobble=I->Wobble;
  p->ramped = ((c1[0]<0.0F)||(c2[0]<0.0F));
  /*
  copy3f(I->WobbleParam,p->wobble_param);*/

  vv=p->v1;
  (*vv++)=(*v1++);
  (*vv++)=(*v1++);
  (*vv++)=(*v1++);
  vv=p->v2;
  (*vv++)=(*v2++);
  (*vv++)=(*v2++);
  (*vv++)=(*v2++);

  I->PrimSize += diff3f(p->v1,p->v2) + 2*r;
  I->PrimSizeCnt++;

  if(I->TTTFlag) {
    transformTTT44f3f(I->TTT,p->v1,p->v1);
    transformTTT44f3f(I->TTT,p->v2,p->v2);
  }

  if(I->Context) {
    RayApplyContextToVertex(I,p->v1);
    RayApplyContextToVertex(I,p->v2);
  }

  vv=p->c1;
  (*vv++)=(*c1++);
  (*vv++)=(*c1++);
  (*vv++)=(*c1++);
  vv=p->c2;
  (*vv++)=(*c2++);
  (*vv++)=(*c2++);
  (*vv++)=(*c2++);
  vv=p->ic;
  {
    float *v;
    vv=p->ic;
    v=I->IntColor;
    (*vv++)=(*v++);
    (*vv++)=(*v++);
    (*vv++)=(*v++);
  }

  I->NPrimitive++;
}
/*========================================================================*/
void RaySausage3fv(CRay *I,float *v1,float *v2,float r,float *c1,float *c2)
{
  CPrimitive *p;

  float *vv;

  VLACacheCheck(I->G,I->Primitive,CPrimitive,I->NPrimitive,0,cCache_ray_primitive);
  p = I->Primitive+I->NPrimitive;

  p->type = cPrimSausage;
  p->r1=r;
  p->trans=I->Trans;
  p->wobble=I->Wobble;
  p->ramped = ((c1[0]<0.0F)||(c2[0]<0.0F));
  /*  
    copy3f(I->WobbleParam,p->wobble_param);*/

  vv=p->v1;
  (*vv++)=(*v1++);
  (*vv++)=(*v1++);
  (*vv++)=(*v1++);
  vv=p->v2;
  (*vv++)=(*v2++);
  (*vv++)=(*v2++);
  (*vv++)=(*v2++);

  I->PrimSize += diff3f(p->v1,p->v2) + 2*r;
  I->PrimSizeCnt++;

  if(I->TTTFlag) {
    transformTTT44f3f(I->TTT,p->v1,p->v1);
    transformTTT44f3f(I->TTT,p->v2,p->v2);
  }

  if(I->Context) {
    RayApplyContextToVertex(I,p->v1);
    RayApplyContextToVertex(I,p->v2);
  }

  vv=p->c1;
  (*vv++)=(*c1++);
  (*vv++)=(*c1++);
  (*vv++)=(*c1++);
  vv=p->c2;
  (*vv++)=(*c2++);
  (*vv++)=(*c2++);
  (*vv++)=(*c2++);
  vv=p->ic;

  {
    float *v=I->IntColor;
    vv=p->ic;
    (*vv++)=(*v++);
    (*vv++)=(*v++);
    (*vv++)=(*v++);
  }

  I->NPrimitive++;
}
void RayEllipsoid3fv(CRay *I,
                     float *v,float r,
                     float *n1,float *n2,float *n3)
{

  CPrimitive *p;

  float *vv;

  VLACacheCheck(I->G,I->Primitive,CPrimitive,I->NPrimitive,0,cCache_ray_primitive);
  p = I->Primitive+I->NPrimitive;
  
  p->type = cPrimEllipsoid;
  p->r1=r; /* maximum extent */
  p->trans=I->Trans;
  p->wobble=I->Wobble;
  p->ramped=(I->CurColor[0]<0.0F);

  I->PrimSize += 2*r;
  I->PrimSizeCnt++;

  vv=p->n0; /* storing lengths of the direction vectors in n0 */

  (*vv++)=length3f(n1);
  (*vv++)=length3f(n2);
  (*vv++)=length3f(n3);

  /* normalize the ellipsoid axes */

  vv=p->n1;
  if(p->n0[0]>R_SMALL8) {
    float factor;
    factor = 1.0F / p->n0[0];
    (*vv++)=(*n1++) * factor;
    (*vv++)=(*n1++) * factor;
    (*vv++)=(*n1++) * factor;
  } else {
    (*vv++)=0.0F;
    (*vv++)=0.0F;
    (*vv++)=0.0F;
  }

  vv=p->n2;
  if(p->n0[1]>R_SMALL8) {
    float factor;
    factor = 1.0F / p->n0[1];
    (*vv++)=(*n2++) * factor;
    (*vv++)=(*n2++) * factor;
    (*vv++)=(*n2++) * factor;
  } else {
    (*vv++)=0.0F;
    (*vv++)=0.0F;
    (*vv++)=0.0F;
  }
  
  vv=p->n3;
  if(p->n0[2]>R_SMALL8) {
    float factor;
    factor = 1.0F / p->n0[2];
    (*vv++)=(*n3++) * factor;
    (*vv++)=(*n3++) * factor;
    (*vv++)=(*n3++) * factor;
  } else {
    (*vv++)=0.0F;
    (*vv++)=0.0F;
    (*vv++)=0.0F;
  }

  vv=p->v1;
  (*vv++)=(*v++);
  (*vv++)=(*v++);
  (*vv++)=(*v++);

  vv=p->c1;
  v=I->CurColor;
  (*vv++)=(*v++);
  (*vv++)=(*v++);
  (*vv++)=(*v++);

  vv=p->ic;
  v=I->IntColor;
  (*vv++)=(*v++);
  (*vv++)=(*v++);
  (*vv++)=(*v++);

  if(I->TTTFlag) {
    transformTTT44f3f(I->TTT,p->v1,p->v1);
    transform_normalTTT44f3f(I->TTT,p->n1,p->n1);
    transform_normalTTT44f3f(I->TTT,p->n2,p->n2);
    transform_normalTTT44f3f(I->TTT,p->n3,p->n3);
  }

  if(I->Context) {
    RayApplyContextToVertex(I,p->v1);
    RayApplyContextToNormal(I,p->n1);
    RayApplyContextToNormal(I,p->n2);
    RayApplyContextToNormal(I,p->n3);
  }

  I->NPrimitive++;

}

/*========================================================================*/
void RayTriangle3fv(CRay *I,
                                      float *v1,float *v2,float *v3,
                                      float *n1,float *n2,float *n3,
                                      float *c1,float *c2,float *c3)
{
  CPrimitive *p;

  float *vv;
  float n0[3],nx[3],s1[3],s2[3],s3[3];
  float l1,l2,l3;

  /*  dump3f(v1," v1");
  dump3f(v2," v2");
  dump3f(v3," v3");
  dump3f(n1," n1");
  dump3f(n2," n2");
  dump3f(n3," n3");
  dump3f(c1," c1");
  dump3f(c2," c2");
  dump3f(c3," c3");*/
  VLACacheCheck(I->G,I->Primitive,CPrimitive,I->NPrimitive,0,cCache_ray_primitive);
  p = I->Primitive+I->NPrimitive;

  p->type = cPrimTriangle;
  p->trans=I->Trans;
  p->tr[0]=I->Trans;
  p->tr[1]=I->Trans;
  p->tr[2]=I->Trans;
  p->wobble=I->Wobble;
  p->ramped = ((c1[0]<0.0F)||(c2[0]<0.0F)||(c3[0]<0.0F));
  /*
    copy3f(I->WobbleParam,p->wobble_param);*/

  /* determine exact triangle normal */
  add3f(n1,n2,nx);
  add3f(n3,nx,nx);
  subtract3f(v1,v2,s1);
  subtract3f(v3,v2,s2);
  subtract3f(v1,v3,s3);
  cross_product3f(s1,s2,n0);
  if((fabs(n0[0])<RAY_SMALL)&&
        (fabs(n0[1])<RAY_SMALL)&&
        (fabs(n0[2])<RAY_SMALL))
       {copy3f(nx,n0);} /* fall-back */
  else if(dot_product3f(n0,nx)<0)
       invert3f(n0);
  normalize3f(n0);

  vv=p->n0;
  (*vv++)=n0[0];
  (*vv++)=n0[1];
  (*vv++)=n0[2];

  /* determine maximum distance from vertex to point */
  l1=(float)length3f(s1);
  l2=(float)length3f(s2);
  l3=(float)length3f(s3);
  if(l2>l1) { if(l3>l2) l1=l3; else l1=l2;  }
  /* store cutoff distance */

  p->r1=l1*0.6F;

  /*  if(l1>20) {
            printf("%8.3f\n",l1);
            printf("%8.3f %8.3f %8.3f\n",s1[0],s1[1],s1[2]);
            printf("%8.3f %8.3f %8.3f\n",s2[0],s2[1],s2[2]);
            printf("%8.3f %8.3f %8.3f\n",s3[0],s3[1],s3[2]);
            }*/

  vv=p->v1;
  (*vv++)=(*v1++);
  (*vv++)=(*v1++);
  (*vv++)=(*v1++);
  vv=p->v2;
  (*vv++)=(*v2++);
  (*vv++)=(*v2++);
  (*vv++)=(*v2++);
  vv=p->v3;
  (*vv++)=(*v3++);
  (*vv++)=(*v3++);
  (*vv++)=(*v3++);

  I->PrimSize += diff3f(p->v1,p->v2) + diff3f(p->v1,p->v3) + diff3f(p->v2,p->v3);
  I->PrimSizeCnt += 3;

  vv=p->c1;
  (*vv++)=(*c1++);
  (*vv++)=(*c1++);
  (*vv++)=(*c1++);
  vv=p->c2;
  (*vv++)=(*c2++);
  (*vv++)=(*c2++);
  (*vv++)=(*c2++);
  vv=p->c3;
  (*vv++)=(*c3++);
  (*vv++)=(*c3++);
  (*vv++)=(*c3++);

  {
    float *v=I->IntColor;
    vv=p->ic;
    (*vv++)=(*v++);
    (*vv++)=(*v++);
    (*vv++)=(*v++);
  }

  vv=p->n1;
  (*vv++)=(*n1++);
  (*vv++)=(*n1++);
  (*vv++)=(*n1++);
  vv=p->n2;
  (*vv++)=(*n2++);
  (*vv++)=(*n2++);
  (*vv++)=(*n2++);
  vv=p->n3;
  (*vv++)=(*n3++);
  (*vv++)=(*n3++);
  (*vv++)=(*n3++);

  if(I->TTTFlag) {
    transformTTT44f3f(I->TTT,p->v1,p->v1);
    transformTTT44f3f(I->TTT,p->v2,p->v2);
    transformTTT44f3f(I->TTT,p->v3,p->v3);
    transform_normalTTT44f3f(I->TTT,p->n0,p->n0);
    transform_normalTTT44f3f(I->TTT,p->n1,p->n1);
    transform_normalTTT44f3f(I->TTT,p->n2,p->n2);
    transform_normalTTT44f3f(I->TTT,p->n3,p->n3);
  }

  if(I->Context) {
    RayApplyContextToVertex(I,p->v1);
    RayApplyContextToVertex(I,p->v2);
    RayApplyContextToVertex(I,p->v3);
    RayApplyContextToNormal(I,p->n0);
    RayApplyContextToNormal(I,p->n1);
    RayApplyContextToNormal(I,p->n2);
    RayApplyContextToNormal(I,p->n3);
  }

  I->NPrimitive++;

}
void RayTriangleTrans3fv(CRay *I,
                         float *v1,float *v2,float *v3,
                         float *n1,float *n2,float *n3,
                         float *c1,float *c2,float *c3,
                         float  t1,float  t2,float  t3)
{
  CPrimitive *p;

  RayTriangle3fv(I,v1,v2,v3,n1,n2,n3,c1,c2,c3);

  p = I->Primitive + I->NPrimitive - 1;
  
  p->tr[0] = t1;
  p->tr[1] = t2;
  p->tr[2] = t3;
  p->trans = (t1+t2+t3)/3.0F;
}

/*========================================================================*/
CRay *RayNew(PyMOLGlobals *G,int antialias)
{
  unsigned int test;
  unsigned char *testPtr;
  int a;

  OOAlloc(I->G,CRay);
  I->G = G;
  test = 0xFF000000;
  testPtr = (unsigned char*)&test;
  I->BigEndian = (*testPtr)&&1;
  I->Trans=0.0F;
  I->Wobble=0;
  I->TTTFlag=false;
  zero3f(I->WobbleParam);
  PRINTFB(I->G,FB_Ray,FB_Blather) 
    " RayNew: BigEndian = %d\n",I->BigEndian
    ENDFB(I->G);

  I->Basis=CacheAlloc(I->G,CBasis,12,0,cCache_ray_basis);
  BasisInit(I->G,I->Basis,0);
  BasisInit(I->G,I->Basis+1,1);
  I->Vert2Prim=VLACacheAlloc(I->G,int,1,0,cCache_ray_vert2prim);
  I->NBasis=2;
  I->Primitive=NULL;
  I->NPrimitive=0;
  I->fColor3fv=RayColor3fv;
  I->fSphere3fv=RaySphere3fv;
  I->fCylinder3fv=RayCylinder3fv;
  I->fCustomCylinder3fv=RayCustomCylinder3fv;
  I->fSausage3fv=RaySausage3fv;
  I->fTriangle3fv=RayTriangle3fv;
  I->fTriangleTrans3fv=RayTriangleTrans3fv;
  I->fCharacter=RayCharacter;
  I->fInteriorColor3fv=RayInteriorColor3fv;
  I->fWobble=RayWobble;
  I->fTransparentf=RayTransparentf;
  I->fEllipsoid3fv=RayEllipsoid3fv;
  I->TTTStackVLA = NULL;
  I->TTTStackDepth = 0;
  I->CheckInterior=false;
  if(antialias<0) antialias = SettingGetGlobal_i(I->G,cSetting_antialias);
  I->Sampling = antialias;
  if(I->Sampling<2) /* always supersample fonts by at least 2X */
    I->Sampling=2;
  for(a=0;a<256;a++) {
    I->Random[a]=(float)((rand()/(1.0+RAND_MAX))-0.5);
  }

  I->Wobble = SettingGet_i(I->G,NULL,NULL,cSetting_ray_texture);
  {
    float *v = SettingGet_3fv(I->G,NULL,NULL,cSetting_ray_texture_settings);
    int color = SettingGetGlobal_color(I->G,cSetting_ray_interior_color);
    copy3f(v,I->WobbleParam);
    v = ColorGet(I->G,color);
    copy3f(v,I->IntColor);
  }

  return(I);
}
/*========================================================================*/
/* BEGIN PROPRIETARY CODE SEGMENT (see disclaimer in "os_proprietary.h") */ 
#ifdef PYMOL_EVAL
#include "RayEvalMessage.h"
#endif
/* END PROPRIETARY CODE SEGMENT */

void RayPrepare(CRay *I,float v0,float v1,float v2,
                float v3,float v4,float v5,
                float fov, float *pos, 
                float *mat,float *rotMat,float aspRat,
                int width, int height, float pixel_scale,int ortho,
                float pixel_ratio,float front_back_ratio,float magnified)
        /*prepare for vertex calls */
{
  int a;
  if(!I->Primitive) 
       I->Primitive=VLACacheAlloc(I->G,CPrimitive,10000,3,cCache_ray_primitive);  
  if(!I->Vert2Prim) 
       I->Vert2Prim=VLACacheAlloc(I->G,int,10000,3,cCache_ray_vert2prim);
  I->Volume[0]=v0;
  I->Volume[1]=v1;
  I->Volume[2]=v2;
  I->Volume[3]=v3;
  I->Volume[4]=v4;
  I->Volume[5]=v5;
  I->Range[0]=I->Volume[1]-I->Volume[0];
  I->Range[1]=I->Volume[3]-I->Volume[2];
  I->Range[2]=I->Volume[5]-I->Volume[4];
  I->AspRatio=aspRat;
  I->Width=width;
  I->Height=height;
  CharacterSetRetention(I->G,true);

  if(mat)  
    for(a=0;a<16;a++)
      I->ModelView[a]=mat[a];
  else {
    for(a=0;a<16;a++)
      I->ModelView[a]=0.0F;
    for(a=0;a<3;a++)
      I->ModelView[a*5]=1.0F;
  }
  if(rotMat)  
    for(a=0;a<16;a++)
      I->Rotation[a]=rotMat[a];
  if(ortho) {
    I->PixelRadius = (((float)I->Range[0])/width)*pixel_scale;
  } else {
    I->PixelRadius = (((float)I->Range[0])/width)*pixel_scale*pixel_ratio;
  }
  I->PixelRatio = pixel_ratio;
  I->Magnified = magnified;
  I->FrontBackRatio = front_back_ratio;
  I->PrimSizeCnt = 0;
  I->PrimSize = 0.0;
  I->Fov = fov;
  copy3f(pos,I->Pos);

/* BEGIN PROPRIETARY CODE SEGMENT (see disclaimer in "os_proprietary.h") */ 
#ifdef PYMOL_EVAL
  RayDrawEvalMessage(I);
#endif
/* END PROPRIETARY CODE SEGMENT */

}
/*========================================================================*/

void RaySetTTT(CRay *I,int flag,float *ttt)
{
  I->TTTFlag=flag;
  if(flag) {
    UtilCopyMem(I->TTT,ttt,sizeof(float)*16);
  }
}

void RayGetTTT(CRay *I,float *ttt)
{
  if(!I->TTTFlag) {
    identity44f(ttt);
  } else {
    copy44f(I->TTT,ttt);
  }
}

/*========================================================================*/
void RayRelease(CRay *I)
{
  int a;

  for(a=0;a<I->NBasis;a++) {
       BasisFinish(&I->Basis[a],a);
  }
  I->NBasis=0;
  VLACacheFreeP(I->G,I->Primitive,0,cCache_ray_primitive,false);
  VLACacheFreeP(I->G,I->Vert2Prim,0,cCache_ray_vert2prim,false);
}
/*========================================================================*/
void RayFree(CRay *I)
{
  RayRelease(I);
  CharacterSetRetention(I->G,false);
  CacheFreeP(I->G,I->Basis,0,cCache_ray_basis,false);
  VLACacheFreeP(I->G,I->Vert2Prim,0,cCache_ray_vert2prim,false);
  VLAFreeP(I->TTTStackVLA);
  OOFreeP(I);
}
/*========================================================================*/
void RayPushTTT(CRay *I)
{
  if(I->TTTFlag) {
    if(!I->TTTStackVLA) {
      I->TTTStackVLA = VLAlloc(float,16);
      copy44f(I->TTT,I->TTTStackVLA);
      I->TTTStackDepth = 1;
    } else {
      float *p;
      VLACheck(I->TTTStackVLA,float,I->TTTStackDepth*16+15);
      p = I->TTTStackVLA + 16*I->TTTStackDepth;
      copy44f(I->TTT,p);
      I->TTTStackDepth++;
    }
  }
}
void RayPopTTT(CRay *I)
{
  if(I->TTTStackDepth>0) {
    float *p;
    I->TTTStackDepth--;
    p = I->TTTStackVLA + 16*I->TTTStackDepth;
    copy44f(p,I->TTT);
    I->TTTFlag = true;
  } else {
    I->TTTFlag = false;
  }
}

void RayApplyMatrix33( unsigned int n, float3 *q, const float m[16],
                          float3 *p )
{
   {
      unsigned int i;
      float m0 = m[0],  m4 = m[4],  m8 = m[8],  m12 = m[12];
      float m1 = m[1],  m5 = m[5],  m9 = m[9],  m13 = m[13];
      float m2 = m[2],  m6 = m[6],  m10 = m[10],  m14 = m[14];
      for (i=0;i<n;i++) {
         float p0 = p[i][0], p1 = p[i][1], p2 = p[i][2];
         q[i][0] = m0 * p0 + m4  * p1 + m8 * p2 + m12;
         q[i][1] = m1 * p0 + m5  * p1 + m9 * p2 + m13;
         q[i][2] = m2 * p0 + m6 * p1 + m10 * p2 + m14;
      }
   }
}

void RayApplyMatrixInverse33( unsigned int n, float3 *q, const float m[16],
                          float3 *p )
{
   {
      unsigned int i;
      float m0 = m[0],  m4 = m[4],  m8 = m[8],  m12 = m[12];
      float m1 = m[1],  m5 = m[5],  m9 = m[9],  m13 = m[13];
      float m2 = m[2],  m6 = m[6],  m10 = m[10],  m14 = m[14];
      for (i=0;i<n;i++) {
         float p0 = p[i][0]-m12, p1 = p[i][1]-m13, p2 = p[i][2]-m14;
         q[i][0] = m0 * p0 + m1  * p1 + m2 * p2;
         q[i][1] = m4 * p0 + m5  * p1 + m6 * p2;
         q[i][2] = m8 * p0 + m9 * p1 + m10 * p2;
      }
   }
}

void RayTransformNormals33( unsigned int n, float3 *q, const float m[16],float3 *p )
{
  unsigned int i;
  float m0 = m[0],  m4 = m[4],  m8 = m[8];
  float m1 = m[1],  m5 = m[5],  m9 = m[9];
  float m2 = m[2],  m6 = m[6],  m10 = m[10];
  for (i=0;i<n;i++) {
    float p0 = p[i][0], p1 = p[i][1], p2 = p[i][2];
    q[i][0] = m0 * p0 + m4  * p1 + m8 * p2;
    q[i][1] = m1 * p0 + m5  * p1 + m9 * p2;
    q[i][2] = m2 * p0 + m6 * p1 + m10 * p2;
  }
  for (i=0;i<n;i++) { /* renormalize - can we do this to the matrix instead? */
    normalize3f(q[i]);
  }
}

void RayTransformInverseNormals33( unsigned int n, float3 *q, const float m[16],float3 *p )
{
  unsigned int i;
  float m0 = m[0],  m4 = m[4],  m8 = m[8];
  float m1 = m[1],  m5 = m[5],  m9 = m[9];
  float m2 = m[2],  m6 = m[6],  m10 = m[10];
  for (i=0;i<n;i++) {
    float p0 = p[i][0], p1 = p[i][1], p2 = p[i][2];
    q[i][0] = m0 * p0 + m1  * p1 + m2 * p2;
    q[i][1] = m4 * p0 + m5  * p1 + m6 * p2;
    q[i][2] = m8 * p0 + m9 * p1 + m10 * p2;
  }
  for (i=0;i<n;i++) { /* renormalize - can we do this to the matrix instead? */
    normalize3f(q[i]);
  }
}




Generated by  Doxygen 1.6.0   Back to index