/* 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) -* -* Z* ------------------------------------------------------------------- */ #ifndef _PYMOL_INLINE #include"os_predef.h" #include"os_std.h" #include"MemoryDebug.h" #include"Base.h" #include"Basis.h" #include"Err.h" #include"Feedback.h" #include"Util.h" #include"MemoryCache.h" #include"Character.h" static const float kR_SMALL4 = 0.0001F; static const float kR_SMALL5 = 0.0001F; #define EPSILON 0.000001 /*========================================================================*/ #ifdef _PYMOL_INLINE __inline__ #endif static int ZLineToSphere(float *base, float *point,float *dir,float radius,float maxial, float *sphere,float *asum,float *pre) { /* Strategy - find an imaginary sphere that lies at the correct point on the line segment, then treat as a sphere reflection */ float intra_p0,intra_p1,intra_p2; float radial,axial,axial_sum,dangle,ab_dangle,axial_perp; float radialsq,tan_acos_dangle; float vradial0,vradial1,vradial2; const float _0 = 0.0f; float point0 = point[0]; float point1 = point[1]; float point2 = point[2]; register float perpAxis0 = pre[0]; /* was cross_product(MinusZ,dir,perpAxis),normalize */ register float perpAxis1 = pre[1]; float intra0 = point0 - base[0]; float intra1 = point1 - base[1]; const float dir0 = dir[0]; const float dir1 = dir[1]; const float dir2 = dir[2]; register float dot, perpDist; /* the perpAxis defines a perp-plane which includes the cyl-axis */ /* get minimum distance between the lines */ perpDist = intra0*perpAxis0 + intra1*perpAxis1; /* was dot_product3f(intra,perpAxis); */ if((perpDist < _0 ? -perpDist : perpDist) > radius) return 0; /*if(fabs(perpDist) > radius) return 0; */ dangle = -dir2; /* was dot(MinusZ,dir) */ ab_dangle = (float)fabs(dangle); if(ab_dangle > (1.0f - kR_SMALL4)) { if(dangle > _0) { sphere[0] = point0; sphere[1] = point1; sphere[2] = point2; } else { sphere[0]=dir0 * maxial + point0; sphere[1]=dir1 * maxial + point1; sphere[2]=dir2 * maxial + point2; } return(1); } if(ab_dangle > kR_SMALL4) tan_acos_dangle = (float)(sqrt1d(1.0-dangle*dangle) / dangle); else tan_acos_dangle = MAXFLOAT; /* now we need to define the triangle in the perp-plane to figure out where the projected line intersection point is */ intra_p2=point2-base[2]; /* first, compute radial distance in the perp-plane between the two starting points */ dot = intra0*perpAxis0 + intra1*perpAxis1; intra_p0=intra0-perpAxis0*dot; intra_p1=intra1-perpAxis1*dot; dot = intra_p0*dir0 + intra_p1*dir1 + intra_p2*dir2; vradial0=intra_p0-dir0*dot; vradial1=intra_p1-dir1*dot; vradial2=intra_p2-dir2*dot; radialsq = ((vradial0*vradial0) + (vradial1*vradial1) + (vradial2*vradial2)); /* now figure out the axial distance along the cyl-line that will give us the point of closest approach */ if(ab_dangle < kR_SMALL4) axial_perp = _0; else axial_perp = (float)(sqrt1f(radialsq)/tan_acos_dangle); axial = (float) sqrt1f( ( (intra_p0*intra_p0) + (intra_p1*intra_p1) + (intra_p2*intra_p2)) - radialsq ); if( (intra_p0*dir0 + intra_p1*dir1 + intra_p2*dir2 ) >= _0 ) axial = axial_perp - axial; else axial = axial_perp + axial; /* now we have to think about where the vector will actually strike the cylinder*/ /* by definition, it must be perpdist away from the perp-plane becuase the perp-plane is parallel to the line, so we can compute the radial component to this point */ radial = radius*radius - perpDist*perpDist; radial = (float)sqrt1f(radial); /* now the trick is figuring out how to adjust the axial distance to get the actual position along the cyl line which will give us a representative sphere */ if(ab_dangle > kR_SMALL4) axial_sum = axial - radial/tan_acos_dangle; else axial_sum = axial; /* printf("radial2 %8.3f \n",radial);*/ if(axial_sum < _0) axial_sum = _0; else if(axial_sum > maxial) axial_sum = maxial; sphere[0] = dir0 * axial_sum + point0; sphere[1] = dir1 * axial_sum + point1; sphere[2] = dir2 * axial_sum + point2; *asum = axial_sum; /* printf("==>%8.3f sphere %8.3f %8.3f %8.3f\n",base[1],sphere[1],axial_perp,axial);*/ return(1); } #ifdef _PYMOL_INLINE __inline__ #endif static int LineToSphere(float *base, float *ray, float *point,float *dir,float radius,float maxial, float *sphere,float *asum) { /* Strategy - find an imaginary sphere that lies at the correct point on the line segment, then treat as a sphere reflection */ register float perpDist,radial,axial,axial_sum,dangle,ab_dangle,axial_perp; register float radialsq,tan_acos_dangle; register float perpAxis0, perpAxis1, perpAxis2; register float intra0, intra1, intra2; register float intra_p0, intra_p1, intra_p2; register float vradial0, vradial1, vradial2; register float dir0 = dir[0], dir1 = dir[1], dir2 = dir[2]; register float ray0 = ray[0], ray1 = ray[1], ray2 = ray[2]; float point0 = point[0], point1 = point[1], point2 = point[2]; register float dot; const float _0 = 0.0F; const float _1 = 1.0F; /* cross_product3f(ray,dir,perpAxis); */ perpAxis0 = ray1*dir2 - ray2*dir1; perpAxis1 = ray2*dir0 - ray0*dir2; perpAxis2 = ray0*dir1 - ray1*dir0; /* subtract3f(point,base,intra); */ intra0 = point0 - base[0]; /* normalize3f(perpAxis) */ { register float len = (float)sqrt1d((perpAxis0 * perpAxis0) + (perpAxis1 * perpAxis1) + (perpAxis2 * perpAxis2)); intra1 = point1 - base[1]; /* subtract3f(point,base,intra); */ if(len>R_SMALL8) { len = _1 / len; intra2 = point2 - base[2]; perpAxis0 *= len; perpAxis1 *= len; perpAxis2 *= len; } else { intra2 = point2 - base[2]; } } /* the perpAxis defines a perp-plane which includes the cyl-axis */ /* get minimum distance between the lines */ /* dot_product3f(intra, perpAxis); */ perpDist = intra0 * perpAxis0; perpDist += intra1 * perpAxis1 + intra2 * perpAxis2 ; /*if(fabs(perpDist) > radius) return 0; */ dangle = ray0 * dir0; if((perpDist < _0 ? -perpDist : perpDist) > radius) return 0; /* dangle = dot_product3f(ray, dir);*/ dangle += ray1 * dir1 + ray2 * dir2; ab_dangle = (float)fabs(dangle); if(ab_dangle > (1.0f - kR_SMALL4)) { if(dangle > _0) { sphere[0] = point0; sphere[1] = point1; sphere[2] = point2; } else { sphere[0]=dir0 * maxial + point0; sphere[1]=dir1 * maxial + point1; sphere[2]=dir2 * maxial + point2; } return(1); } if(ab_dangle > kR_SMALL4) tan_acos_dangle = (float)(sqrt1d(1.0-dangle*dangle) / dangle); else tan_acos_dangle = MAXFLOAT; /* now we need to define the triangle in the perp-plane to figure out where the projected line intersection point is */ /* first, compute radial distance in the perp-plane between the two starting points */ /* dot = dot_product3f(intra,perpAxis); */ dot = intra0 * perpAxis0 + intra1 * perpAxis1 + intra2 * perpAxis2 ; intra_p0 = intra0-perpAxis0*dot; intra_p1 = intra1-perpAxis1*dot; intra_p2 = intra2-perpAxis2*dot; /* dot = dot_product3f(intra_p, dir); */ dot = intra_p0 * dir0 + intra_p1 * dir1 + intra_p2 * dir2; vradial0=intra_p0-dir0*dot; vradial1=intra_p1-dir1*dot; vradial2=intra_p2-dir2*dot; radialsq = ((vradial0*vradial0) + (vradial1*vradial1) + (vradial2*vradial2)); /* now figure out the axial distance along the cyl-line that will give us the point of closest approach */ if(ab_dangle < kR_SMALL4) axial_perp = _0; else axial_perp = (float)(sqrt1f(radialsq)/tan_acos_dangle); axial = (float) sqrt1f( ( (intra_p0*intra_p0) + (intra_p1*intra_p1) + (intra_p2*intra_p2)) - radialsq ); if( (intra_p0*dir0 + intra_p1*dir1 + intra_p2*dir2 ) >= _0 ) axial = axial_perp - axial; else axial = axial_perp + axial; /* now we have to think about where the vector will actually strike the cylinder*/ /* by definition, it must be perpdist away from the perp-plane becuase the perp-plane is parallel to the line, so we can compute the radial component to this point */ radial = radius*radius - perpDist*perpDist; radial = (float)sqrt1f(radial); /* now the trick is figuring out how to adjust the axial distance to get the actual position along the cyl line which will give us a representative sphere */ if(ab_dangle > kR_SMALL4) axial_sum = axial - radial/tan_acos_dangle; else axial_sum = axial; /* printf("radial2 %8.3f \n",radial);*/ if(axial_sum < _0) axial_sum = _0; else if(axial_sum > maxial) axial_sum = maxial; sphere[0] = dir0 * axial_sum + point0; sphere[1] = dir1 * axial_sum + point1; sphere[2] = dir2 * axial_sum + point2; *asum = axial_sum; /* printf("==>%8.3f sphere %8.3f %8.3f %8.3f\n",base[1],sphere[1],axial_perp,axial);*/ return(1); } #ifdef _PYMOL_INLINE __inline__ #endif static int FrontToInteriorSphere(float *front, float *point, float *dir, float radius, float radius2, float maxial) { float intra_p[3]; float axial; float intra[3],axis[3]; float sphere[3]; subtract3f(point,front,intra); remove_component3f(intra,dir,intra_p); add3f(front,intra_p,intra_p); subtract3f(point,intra_p,axis); axial = -dot_product3f(axis,dir); if(axial<0.0F) axial=0.0F; else if(axial>maxial) axial = maxial; sphere[0]=axial*dir[0]+point[0]; sphere[1]=axial*dir[1]+point[1]; sphere[2]=axial*dir[2]+point[2]; return (diffsq3f(sphere,front)<=radius2); } /*========================================================================*/ #ifdef _PYMOL_INLINE __inline__ #endif static int ZLineToSphereCapped(float *base,float *point, float *dir,float radius,float maxial, float *sphere,float *asum,int cap1, int cap2,float *pre) { /* Strategy - find an imaginary sphere that lies at the correct point on the line segment, then treat as a sphere reflection */ float perpAxis[3],intra_p[3]; float perpDist,radial,axial,axial_sum,dangle,ab_dangle,axial_perp; float radialsq,tan_acos_dangle; float len_proj; float intra[3],vradial[3]; float diff[3],fpoint[3]; float proj[3]; perpAxis[0] = pre[0]; /* was cross_product(MinusZ,dir,perpAxis),normalize */ perpAxis[1] = pre[1]; /* the perpAxis defines a perp-plane which includes the cyl-axis */ intra[0]=point[0]-base[0]; intra[1]=point[1]-base[1]; /* get minimum distance between the lines */ perpDist = intra[0]*perpAxis[0] + intra[1]*perpAxis[1]; /* was dot_product3f(intra,perpAxis); */ if(fabs(perpDist)>radius) { return(0); } perpAxis[2] = 0.0; intra[2]=point[2]-base[2]; dangle = -dir[2]; /* was dot(MinusZ,dir) */ ab_dangle= (float)fabs(dangle); if(ab_dangle>(1-kR_SMALL4)) /* vector inline with light ray... */ { vradial[0]=point[0]-base[0]; vradial[1]=point[1]-base[1]; vradial[2]=0.0F; radial = (float)length3f(vradial); if(radial>radius) return 0; if(dangle>0.0) { switch(cap1) { case cCylCapFlat: sphere[0]=base[0]; sphere[1]=base[1]; sphere[2]=point[2]-radius; break; case cCylCapRound: sphere[0]=point[0]; sphere[1]=point[1]; sphere[2]=point[2]; } return(1); } else { switch(cap1) { case cCylCapFlat: sphere[0]=base[0]; sphere[1]=base[1]; sphere[2]=dir[2]*maxial+point[2]-radius; break; case cCylCapRound: sphere[0]=dir[0]*maxial+point[0]; sphere[1]=dir[1]*maxial+point[1]; sphere[2]=dir[2]*maxial+point[2]; } return(1); } } /*tan_acos_dangle = tan(acos(dangle));*/ tan_acos_dangle = (float)sqrt1f(1-dangle*dangle)/dangle; /* printf("perpDist %8.3f\n",perpDist); printf("dir %8.3f %8.3f %8.3f\n",dir[0],dir[1],dir[2]); printf("base %8.3f %8.3f %8.3f\n",base[0],base[1],base[2]); printf("point %8.3f %8.3f %8.3f\n",point[0],point[1],point[2]); printf("intra %8.3f %8.3f %8.3f\n",intra[0],intra[1],intra[2]); printf("perpAxis %8.3f %8.3f %8.3f\n",perpAxis[0],perpAxis[1],perpAxis[2]); */ /* now we need to define the triangle in the perp-plane to figure out where the projected line intersection point is */ /* first, compute radial distance in the perp-plane between the two starting points */ remove_component3f(intra,perpAxis,intra_p); remove_component3f(intra_p,dir,vradial); radialsq = lengthsq3f(vradial); /* now figure out the axial distance along the cyl-line that will give us the point of closest approach */ if(ab_dangle<kR_SMALL4) axial_perp=0; else axial_perp = (float)sqrt1f(radialsq)/tan_acos_dangle; axial = (float)lengthsq3f(intra_p)-radialsq; axial = (float)sqrt1f(axial); /* printf("radial %8.3f\n",radial); printf("vradial %8.3f %8.3f %8.3f\n",vradial[0],vradial[1],vradial[2]); printf("radial %8.3f\n",radial); printf("dangle %8.3f \n",dangle); printf("axial_perp %8.3f \n",axial_perp); printf("axial1 %8.3f \n",axial); printf("%8.3f\n",dot_product3f(intra_p,dir)); */ if(dot_product3f(intra_p,dir)>=0.0) axial = axial_perp - axial; else axial = axial_perp + axial; /* printf("axial2 %8.3f\n",axial); */ /* now we have to think about where the vector will actually strike the cylinder*/ /* by definition, it must be perpdist away from the perp-plane becuase the perp-plane is parallel to the line, so we can compute the radial component to this point */ radial = radius*radius-perpDist*perpDist; radial = (float)sqrt1f(radial); /* now the trick is figuring out how to adjust the axial distance to get the actual position along the cyl line which will give us a representative sphere */ if(ab_dangle > kR_SMALL4) axial_sum = axial - radial/tan_acos_dangle; else axial_sum = axial; /* printf("ab_dangle %8.3f \n",ab_dangle); printf("axial_sum %8.3f \n",axial_sum); */ if(axial_sum<0) { switch(cap1) { case cCylCapFlat: subtract3f(point,base,diff); project3f(diff,dir,proj); len_proj = (float)length3f(proj); dangle = -proj[2]/len_proj; if(fabs(dangle)<kR_SMALL4) return 0; sphere[0]=base[0]; sphere[1]=base[1]; sphere[2]=base[2]-len_proj/dangle; if(diff3f(sphere,point)>radius) return 0; sphere[0]+=dir[0]*radius; sphere[1]+=dir[1]*radius; sphere[2]+=dir[2]*radius; *asum=0; break; case cCylCapRound: axial_sum=0; sphere[0]=point[0]; sphere[1]=point[1]; sphere[2]=point[2]; *asum = axial_sum; break; case cCylCapNone: default: return 0; break; } } else if(axial_sum>maxial) { switch(cap2) { case cCylCapFlat: scale3f(dir,maxial,fpoint); add3f(fpoint,point,fpoint); subtract3f(fpoint,base,diff); project3f(diff,dir,proj); len_proj = (float)length3f(proj); dangle = -proj[2]/len_proj; if(fabs(dangle)<kR_SMALL4) return 0; sphere[0]=base[0]; sphere[1]=base[1]; sphere[2]=base[2]-len_proj/dangle; if(diff3f(sphere,fpoint)>radius) return 0; sphere[0]-=dir[0]*radius; sphere[1]-=dir[1]*radius; sphere[2]-=dir[2]*radius; *asum=maxial; break; case cCylCapRound: axial_sum=maxial; sphere[0]=dir[0]*axial_sum+point[0]; sphere[1]=dir[1]*axial_sum+point[1]; sphere[2]=dir[2]*axial_sum+point[2]; *asum = axial_sum; break; case cCylCapNone: default: return 0; break; } } else { sphere[0]=dir[0]*axial_sum+point[0]; sphere[1]=dir[1]*axial_sum+point[1]; sphere[2]=dir[2]*axial_sum+point[2]; *asum = axial_sum; /* printf("==>%8.3f sphere %8.3f %8.3f %8.3f\n",base[1],sphere[1],axial_perp,axial);*/ } return(1); } #ifdef _PYMOL_INLINE __inline__ #endif static int LineToSphereCapped(float *base, float *ray, float *point,float *dir,float radius,float maxial, float *sphere,float *asum,int cap1,int cap2) { /* Strategy - find an imaginary sphere that lies at the correct point on the line segment, then treat as a sphere reflection */ float perpAxis[3],intra_p[3]; float perpDist,radial,axial,axial_sum,dangle,ab_dangle,axial_perp; float radialsq,tan_acos_dangle; float len_proj; float intra[3],vradial[3]; float diff[3],fpoint[3]; float proj[3]; subtract3f(point,base,intra); cross_product3f(ray,dir,perpAxis); normalize3f(perpAxis); /* the perpAxis defines a perp-plane which includes the cyl-axis */ /* get minimum distance between the lines */ perpDist = dot_product3f(intra, perpAxis); /* was dot_product3f(intra,perpAxis); */ if(fabs(perpDist)>radius) { return(0); } dangle = dot_product3f(ray, dir); ab_dangle = (float)fabs(dangle); if(ab_dangle>(1-kR_SMALL4)) /* vector inline with light ray... */ { vradial[0]=point[0]-base[0]; vradial[1]=point[1]-base[1]; vradial[2]=point[2]-base[2]; radial = (float)length3f(vradial); if(radial>radius) return 0; if(dangle>0.0) { switch(cap1) { case cCylCapFlat: sphere[0]=dir[0]*radius+point[0]; sphere[1]=dir[1]*radius+point[1]; sphere[2]=dir[2]*radius+point[2]; break; case cCylCapRound: sphere[0]=point[0]; sphere[1]=point[1]; sphere[2]=point[2]; } return(1); } else { switch(cap1) { case cCylCapFlat: maxial-=radius; break; } sphere[0]=dir[0]*maxial+point[0]; sphere[1]=dir[1]*maxial+point[1]; sphere[2]=dir[2]*maxial+point[2]; return(1); } } /*tan_acos_dangle = tan(acos(dangle));*/ tan_acos_dangle = (float)sqrt1f(1-dangle*dangle)/dangle; /* printf("perpDist %8.3f\n",perpDist); printf("dir %8.3f %8.3f %8.3f\n",dir[0],dir[1],dir[2]); printf("base %8.3f %8.3f %8.3f\n",base[0],base[1],base[2]); printf("point %8.3f %8.3f %8.3f\n",point[0],point[1],point[2]); printf("intra %8.3f %8.3f %8.3f\n",intra[0],intra[1],intra[2]); printf("perpAxis %8.3f %8.3f %8.3f\n",perpAxis[0],perpAxis[1],perpAxis[2]); */ /* now we need to define the triangle in the perp-plane to figure out where the projected line intersection point is */ /* first, compute radial distance in the perp-plane between the two starting points */ remove_component3f(intra,perpAxis,intra_p); remove_component3f(intra_p,dir,vradial); radialsq = lengthsq3f(vradial); /* now figure out the axial distance along the cyl-line that will give us the point of closest approach */ if(ab_dangle<kR_SMALL4) axial_perp=0; else axial_perp = (float)sqrt1f(radialsq)/tan_acos_dangle; axial = (float)lengthsq3f(intra_p)-radialsq; axial = (float)sqrt1f(axial); /* printf("radial %8.3f\n",radial); printf("vradial %8.3f %8.3f %8.3f\n",vradial[0],vradial[1],vradial[2]); printf("radial %8.3f\n",radial); printf("dangle %8.3f \n",dangle); printf("axial_perp %8.3f \n",axial_perp); printf("axial1 %8.3f \n",axial); printf("%8.3f\n",dot_product3f(intra_p,dir)); */ if(dot_product3f(intra_p,dir)>=0.0) axial = axial_perp - axial; else axial = axial_perp + axial; /* printf("axial2 %8.3f\n",axial); */ /* now we have to think about where the vector will actually strike the cylinder*/ /* by definition, it must be perpdist away from the perp-plane becuase the perp-plane is parallel to the line, so we can compute the radial component to this point */ radial = radius*radius-perpDist*perpDist; radial = (float)sqrt1f(radial); /* now the trick is figuring out how to adjust the axial distance to get the actual position along the cyl line which will give us a representative sphere */ if(ab_dangle > kR_SMALL4) axial_sum = axial - radial/tan_acos_dangle; else axial_sum = axial; /* printf("ab_dangle %8.3f \n",ab_dangle); printf("axial_sum %8.3f \n",axial_sum); */ if(axial_sum<0) { switch(cap1) { case cCylCapFlat: subtract3f(point,base,diff); project3f(diff,dir,proj); len_proj = (float)length3f(proj); dangle = dot_product3f(proj,ray)/len_proj; if(fabs(dangle)<kR_SMALL4) return 0; len_proj /= dangle; sphere[0]=base[0]+ray[0]*len_proj; sphere[1]=base[1]+ray[1]*len_proj; sphere[2]=base[2]+ray[2]*len_proj; if(diff3f(sphere,point)>radius) return 0; sphere[0]+=dir[0]*radius; sphere[1]+=dir[1]*radius; sphere[2]+=dir[2]*radius; *asum=0; break; case cCylCapRound: axial_sum=0; sphere[0]=point[0]; sphere[1]=point[1]; sphere[2]=point[2]; *asum = axial_sum; break; case cCylCapNone: default: return 0; break; } } else if(axial_sum>maxial) { switch(cap2) { case cCylCapFlat: scale3f(dir,maxial,fpoint); add3f(fpoint,point,fpoint); subtract3f(fpoint,base,diff); project3f(diff,dir,proj); len_proj = (float)length3f(proj); dangle = dot_product3f(proj,ray)/len_proj; if(fabs(dangle)<kR_SMALL4) return 0; len_proj /= dangle; sphere[0]=base[0]+ray[0]*len_proj; sphere[1]=base[1]+ray[1]*len_proj; sphere[2]=base[2]+ray[2]*len_proj; if(diff3f(sphere,fpoint)>radius) return 0; sphere[0]-=dir[0]*radius; sphere[1]-=dir[1]*radius; sphere[2]-=dir[2]*radius; *asum=maxial; break; case cCylCapRound: axial_sum=maxial; sphere[0]=dir[0]*axial_sum+point[0]; sphere[1]=dir[1]*axial_sum+point[1]; sphere[2]=dir[2]*axial_sum+point[2]; *asum = axial_sum; break; case cCylCapNone: default: return 0; break; } } else { sphere[0]=dir[0]*axial_sum+point[0]; sphere[1]=dir[1]*axial_sum+point[1]; sphere[2]=dir[2]*axial_sum+point[2]; *asum = axial_sum; /* printf("==>%8.3f sphere %8.3f %8.3f %8.3f\n",base[1],sphere[1],axial_perp,axial);*/ } return(1); } static int ConeLineToSphereCapped(float *base, float *ray, float *point,float *dir,float radius, float small_radius, float maxial, float *sphere, float *asum, float *sph_rad, float *sph_rad_sq, int cap1,int cap2) { /* Strategy - find an imaginary sphere that lies at the correct point on the line segment, then treat as a sphere reflection */ float axial_sum, dangle, ab_dangle; float len_proj; float diff[3],fpoint[3]; float proj[3]; { float perp_axis[3]; float intra[3]; float perp_dist; subtract3f(point,base,intra); cross_product3f(ray,dir,perp_axis); normalize3f(perp_axis); /* the perp_axis defines a perp-plane which includes the cyl-axis */ /* get minimum distance between the lines */ perp_dist = fabs(dot_product3f(intra, perp_axis)); if(perp_dist>radius) { /* the infinite ray and the cone direction lines don't pass close enough to intersect within the bounding cylinder */ return 0; } } dangle = dot_product3f(ray, dir); ab_dangle = (float)fabs(dangle); /* set up the cone */ { double spread = (radius - small_radius) / maxial; float orig_axial_len = radius / spread; float base2orig_radial[3]; float base2orig_normal[3]; float base2orig_radial_len, base2orig_radial_len_sq; float base2orig_len, base2orig_len_sq; float base2orig_axial_len, base2orig_spread; float orig[3]; float base2orig[3]; float near_p[3]; int base_inside_cone = false; float shift1, shift2; /* this is what we are solving for -- the distance along the cone axis to the point of intersection */ scale3f(dir, orig_axial_len, orig); add3f(point, orig, orig); subtract3f(orig, base ,base2orig); remove_component3f(base2orig, dir, base2orig_radial); base2orig_radial_len_sq = lengthsq3f(base2orig_radial); base2orig_len_sq = lengthsq3f(base2orig); base2orig_axial_len = sqrt1f(base2orig_len_sq - base2orig_radial_len_sq); base2orig_radial_len = sqrt1f(base2orig_radial_len_sq); base2orig_len = sqrt1f(base2orig_len_sq); base2orig_spread = base2orig_radial_len / base2orig_axial_len; base_inside_cone = (base2orig_spread < spread); normalize23f(base2orig, base2orig_normal); if(ab_dangle > kR_SMALL4) { float ray_extend = base2orig_axial_len / dangle; if(dot_product3f(base2orig_normal,dir)<0.0) ray_extend = -ray_extend; scale3f(ray, ray_extend, near_p); add3f(base, near_p, near_p); #if 0 printf("base2orig len %8.3f extend %8.3f dangle %8.3f %8.3f ", base2orig_radial_len,ray_extend,dangle, dot_product3f(base2orig_normal, dir)); { float check[3]; subtract3f(near_p,orig,check); printf("check: %8.7f\n",dot_product3f(check,dir)); } #endif /* Now we punt entirely and throw the solution of this quadratic relationship over to Mathematica. Surely this calculation could be significantly optimized... */ { double partA, partB, partC; double dir0 = dir[0], dir1 = dir[1], dir2 = dir[2]; double ray0 = ray[0], ray1 = ray[1], ray2 = ray[2]; double dir0Sq = dir0*dir0, dir1Sq = dir1*dir1, dir2Sq = dir2*dir2; double ray0Sq = ray0*ray0, ray1Sq = ray1*ray1, ray2Sq = ray2*ray2; double cone0 = orig[0], cone1 = orig[1], cone2 = orig[2]; double near0 = near_p[0], near1 = near_p[1], near2 = near_p[2]; double cone0Sq = cone0*cone0, cone1Sq = cone1*cone1, cone2Sq = cone2*cone2; double near0Sq = near0*near0, near1Sq = near1*near1, near2Sq = near2*near2; double dAngle = ray0 * dir0 + ray1 * dir1 + ray2 * dir2; double spreadSq = spread * spread; double dAngleSq = dAngle * dAngle; partB = dAngleSq*(4*pow(cone0*dAngle*dir0 + cone1*dAngle*dir1 + cone2*dAngle*dir2 - dAngle*dir0*near0 - dAngle*dir1*near1 - dAngle*dir2*near2 - cone0*ray0 + near0*ray0 - cone1*ray1 + near1*ray1 - cone2*ray2 + near2*ray2,2.0) - 4*(cone0Sq + cone1Sq + cone2Sq - 2*cone0*near0 + near0Sq - 2*cone1*near1 + near1Sq - 2*cone2*near2 + near2Sq)* (ray0Sq + ray1Sq - 2*dAngle*(dir0*ray0 + dir1*ray1 + dir2*ray2) + ray2Sq + dAngleSq*(dir0Sq + dir1Sq + dir2Sq - spreadSq))); if(partB<0.0) { /* negative? then there are NO real solutions */ return 0; } else { double partBroot = sqrt(partB); partA = -(cone0*dAngleSq*dir0) - cone1*dAngleSq*dir1 - cone2*dAngleSq*dir2 + dAngleSq*dir0*near0 + dAngleSq*dir1*near1 + dAngleSq*dir2*near2 + cone0*dAngle*ray0 - dAngle*near0*ray0 + cone1*dAngle*ray1 - dAngle*near1*ray1 + cone2*dAngle*ray2 - dAngle*near2*ray2; partC = (ray0Sq + ray1Sq - 2*dAngle*(dir0*ray0 + dir1*ray1 + dir2*ray2) + ray2Sq + dAngleSq*(dir0Sq + dir1Sq + dir2Sq - spreadSq)); shift1 = (float)((partA + partBroot*0.5) / partC); shift2 = (float)((partA - partBroot*0.5) / partC); } #if 0 /* prove that the solutions are valid */ { float near1[3],cone1[3]; float near2[3],cone2[3]; float radius1 = fabs(spread * shift1); float radius2 = fabs(spread * shift2); scale3f(ray,shift1/dangle,near1); add3f(near_p,near1,near1); scale3f(dir,shift1,cone1); add3f(orig,cone1,cone1); scale3f(ray,shift2/dangle,near2); add3f(near_p,near2,near2); scale3f(dir,shift2,cone2); add3f(orig,cone2,cone2); { float diff1 = diff3f(near1,cone1); float diff2 = diff3f(near2,cone2); if(fabs(diff1 - radius1)>0.01) { printf("error #1: %10.8f != %10.8f for ", diff1, radius1); printf("%8.3f %8.3f\n",shift1,shift2); } if(fabs(diff2 - radius2)>0.01) { printf("error #2: %10.8f != %10.8f for ", diff2, radius2); printf("%8.3f %8.3f\n",shift1,shift2); } } } #endif } { float axial_sum1 = orig_axial_len + shift1; float axial_sum2 = orig_axial_len + shift2; if(dangle>0.0F) { /* cone is narrowing in parallel with ray */ if(shift1<shift2) { axial_sum = axial_sum1; } else { axial_sum = axial_sum2; } if((axial_sum<0.0F) || (base_inside_cone && (axial_sum<orig_axial_len))) { switch(cap1) { case cCylCapFlat: subtract3f(point,base,diff); project3f(diff,dir,proj); len_proj = (float)length3f(proj); dangle = dot_product3f(proj,ray)/len_proj; if(fabs(dangle)<kR_SMALL5) return 0; len_proj /= dangle; sphere[0]=base[0]+ray[0]*len_proj; sphere[1]=base[1]+ray[1]*len_proj; sphere[2]=base[2]+ray[2]*len_proj; if(diff3f(sphere,point)>radius) return 0; sphere[0]+=dir[0]*radius; sphere[1]+=dir[1]*radius; sphere[2]+=dir[2]*radius; *sph_rad = radius; *sph_rad_sq = radius * radius; *asum=0; return 1; break; case cCylCapNone: default: return 0; break; } } else if(axial_sum>maxial) { return 0; } } else { /* cone is narrowing against ray */ if(shift1<shift2) { axial_sum = axial_sum2; if(axial_sum>orig_axial_len) { axial_sum = axial_sum1; } } else { axial_sum = axial_sum1; if(axial_sum>orig_axial_len) { axial_sum = axial_sum2; } } if(axial_sum<0.0F) { return 0; } else if(axial_sum>=maxial) { switch(cap2) { case cCylCapFlat: scale3f(dir,maxial,fpoint); add3f(fpoint,point,fpoint); subtract3f(fpoint,base,diff); project3f(diff,dir,proj); len_proj = (float)length3f(proj); dangle = dot_product3f(proj,ray)/len_proj; if(fabs(dangle)<kR_SMALL5) return 0; len_proj /= dangle; sphere[0]=base[0]+ray[0]*len_proj; sphere[1]=base[1]+ray[1]*len_proj; sphere[2]=base[2]+ray[2]*len_proj; if(diff3f(sphere,fpoint)>small_radius) /* need to handle this case */ return 0; sphere[0]-=dir[0]*small_radius; sphere[1]-=dir[1]*small_radius; sphere[2]-=dir[2]*small_radius; *sph_rad = small_radius; *sph_rad_sq = small_radius * small_radius; *asum=maxial; return 1; break; case cCylCapNone: default: return 0; break; } } } } } else { axial_sum = orig_axial_len - base2orig_axial_len; if((axial_sum<0.0F)||(axial_sum>maxial)) return 0; } /* normal hit in mid-section of the cone */ { float radius_at_hit = radius - spread * axial_sum; float adjustment = radius_at_hit * spread; *asum = axial_sum; /* color blend based on actual hit location */ axial_sum -= adjustment; /* printf(" %8.3f ",axial_sum);*/ sphere[0] = dir[0] * axial_sum + point[0]; sphere[1] = dir[1] * axial_sum + point[1]; sphere[2] = dir[2] * axial_sum + point[2]; /* and provide virtual sphere radius info */ *sph_rad_sq = (float)((radius_at_hit * radius_at_hit) + (adjustment * adjustment)); *sph_rad = (float)sqrt(*sph_rad_sq); /* printf("%8.3f %8.3f ",*sph_rad, tan_acos_dangle);*/ /* dump3f(sphere,"sph");*/ return 1; } } return 0; } #ifdef _PYMOL_INLINE __inline__ #endif static int FrontToInteriorSphereCapped(float *front, float *point, float *dir, float radius, float radius2, float maxial, int cap1, int cap2) { float intra_p[3]; float axial; float intra[3],axis[3]; float sphere[3]; subtract3f(point,front,intra); remove_component3f(intra,dir,intra_p); add3f(front,intra_p,intra_p); subtract3f(point,intra_p,axis); axial = -dot_product3f(axis,dir); if(axial<0.0F) return 0; else if(axial>maxial) return 0; sphere[0]=axial*dir[0]+point[0]; sphere[1]=axial*dir[1]+point[1]; sphere[2]=axial*dir[2]+point[2]; return (diffsq3f(sphere,front)<radius2); } /*========================================================================*/ #ifdef _PYMOL_INLINE __inline__ #endif static float ZLineClipPoint(float *base,float *point,float *alongNormalSq,float cutoff) { float hyp0, hyp1, hyp2; float result = MAXFLOAT; /* this routine determines whether or not a vector starting at "base" heading in the direction "normal" intersects a sphere located at "point". It returns how far along the vector the intersection with the plane is or MAXFLOAT if there isn't a relevant intersection NOTE: this routine has been optimized for normals along Z Optimizes-out vectors that are more than "cutoff" from "point" in x,y plane */ hyp0 = point[0] - base[0]; if(fabs(hyp0) > cutoff) return result; hyp1 = point[1] - base[1]; if(fabs(hyp1) > cutoff) return result; hyp2 = point[2] - base[2]; if(hyp2 < 0.0) { (*alongNormalSq) = (hyp2 * hyp2); result = (hyp0 * hyp0) + (hyp1 * hyp1); } return result; } #ifdef _PYMOL_INLINE __inline__ #endif static float ZLineClipPointNoZCheck(float *base,float *point,float *alongNormalSq,float cutoff) { float hyp0, hyp1, hyp2; float result = MAXFLOAT; /* this routine determines whether or not a vector starting at "base" heading in the direction "normal" intersects a sphere located at "point". It returns how far along the vector the intersection with the plane is or MAXFLOAT if there isn't a relevant intersection NOTE: this routine has been optimized for normals along Z Optimizes-out vectors that are more than "cutoff" from "point" in x,y plane */ hyp0 = point[0] - base[0]; if(fabs(hyp0) > cutoff) return result; hyp1 = point[1] - base[1]; if(fabs(hyp1) > cutoff) return result; hyp2 = point[2] - base[2]; (*alongNormalSq) = (hyp2 * hyp2); result = (hyp0 * hyp0) + (hyp1 * hyp1); return result; } #ifdef _PYMOL_INLINE __inline__ #endif static int LineClipPoint(float *base,float *ray, float *point, float *dist, float cutoff, float cutoff2) { register float hyp0, hyp1, hyp2; register float opp0, opp1, opp2; register float adj0, adj1, adj2; register float ray0, ray1, ray2; register float proj; register double dcutoff = (double)cutoff; register float opp_len_sq; /* this routine determines whether or not a vector starting at "base" heading in the direction "ray" intersects a sphere located at "point". It returns how far along the vector the intersection with the plane is or MAXFLOAT if there isn't a relevant intersection NOTE: this routine has been optimized for normals along Z Optimizes-out vectors that are more than "cutoff" from "point" */ /* compute the hypo */ hyp2 = point[2] - base[2]; hyp1 = point[1] - base[1]; hyp0 = point[0] - base[0]; ray0 = ray[0]; ray1 = ray[1]; ray2 = ray[2]; /* compute the adjacent edge (dot-projection) */ proj = (ray0 * hyp0) + (ray1 * hyp1) + (ray2 * hyp2); adj0 = ray0 * proj; adj1 = ray1 * proj; opp0 = hyp0 - adj0; adj2 = ray2 * proj; if(fabs(opp0) > dcutoff) return 0; opp1 = hyp1 - adj1; opp2 = hyp2 - adj2; if(fabs(opp1) > dcutoff) return 0; if(fabs(opp2) > dcutoff) return 0; opp_len_sq = (opp0 * opp0) + (opp1 * opp1) + (opp2 * opp2); if(opp_len_sq<=cutoff2) { *dist = proj - (float)sqrt1f(cutoff2-opp_len_sq); return 1; } return 0; } static int LineClipEllipsoidPoint(float *base,float *ray, float *point, float *dist, float cutoff, float cutoff2, float *scale, float *n1, float *n2, float *n3) { float point_to_base[3]; float scaled_base[3]; float scaled_ray[3]; float d1,d2,d3,s1,s2,s3; float comp1[3], comp2[3], comp3[3]; subtract3f(base, point, point_to_base); /* project difference vector onto ellipsoid axes */ d1 = dot_product3f(point_to_base, n1); d2 = dot_product3f(point_to_base, n2); d3 = dot_product3f(point_to_base, n3); s1 = d1 / scale[0]; s2 = d2 / scale[1]; s3 = d3 / scale[2]; scale3f(n1, s1, comp1); scale3f(n2, s2, comp2); scale3f(n3, s3, comp3); copy3f(point, scaled_base); add3f(comp1, scaled_base, scaled_base); add3f(comp2, scaled_base, scaled_base); add3f(comp3, scaled_base, scaled_base); d1 = dot_product3f(ray, n1); d2 = dot_product3f(ray, n2); d3 = dot_product3f(ray, n3); s1 = d1 / scale[0]; s2 = d2 / scale[1]; s3 = d3 / scale[2]; scale3f(n1, s1, comp1); scale3f(n2, s2, comp2); scale3f(n3, s3, comp3); copy3f(comp1, scaled_ray); add3f(comp2, scaled_ray, scaled_ray); add3f(comp3, scaled_ray, scaled_ray); d1 = length3f(scaled_ray); normalize3f(scaled_ray); /* dump3f(ray,"ray"); dump3f(scaled_ray,"scaled_ray"); */ { register float hyp0, hyp1, hyp2; register float opp0, opp1, opp2; register float adj0, adj1, adj2; register float ray0, ray1, ray2; register float proj; register double dcutoff = (double)cutoff; register float opp_len_sq; /* this routine determines whether or not a vector starting at "base" heading in the direction "ray" intersects a sphere located at "point". It returns how far along the vector the intersection with the plane is or MAXFLOAT if there isn't a relevant intersection */ /* compute the hypo */ hyp2 = point[2] - scaled_base[2]; hyp1 = point[1] - scaled_base[1]; hyp0 = point[0] - scaled_base[0]; ray0 = scaled_ray[0]; ray1 = scaled_ray[1]; ray2 = scaled_ray[2]; /* compute the adjacent edge (dot-projection) */ proj = (ray0 * hyp0) + (ray1 * hyp1) + (ray2 * hyp2); adj0 = ray0 * proj; adj1 = ray1 * proj; opp0 = hyp0 - adj0; adj2 = ray2 * proj; if(fabs(opp0) > dcutoff) return 0; opp1 = hyp1 - adj1; opp2 = hyp2 - adj2; if(fabs(opp1) > dcutoff) return 0; if(fabs(opp2) > dcutoff) return 0; opp_len_sq = (opp0 * opp0) + (opp1 * opp1) + (opp2 * opp2); if(opp_len_sq<=cutoff2) { /* line hits the virtual sphere */ float scaled_dist = proj - (float)sqrt1f(cutoff2-opp_len_sq); *(dist) = scaled_dist / d1; return 1; } return 0; } } /*========================================================================*/ void BasisSetupMatrix(CBasis *I) { float oldZ[3] = { 0.0,0.0,1.0 }; float newY[3]; float dotgle,angle; cross_product3f(oldZ,I->LightNormal,newY); dotgle=dot_product3f(oldZ,I->LightNormal); if((1.0-fabs(dotgle))<kR_SMALL4) { dotgle=(float)(dotgle/fabs(dotgle)); newY[0]=0.0; newY[1]=1.0; newY[2]=0.0; } normalize3f(newY); angle=(float)(-acos(dotgle)); /* now all we gotta do is effect a rotation about the new Y axis to line up new Z with Z */ rotation_to_matrix33f(newY,angle,I->Matrix); /* printf("%8.3f %8.3f %8.3f %8.3f\n",angle*180.0/cPI,newY[0],newY[1],newY[2]); matrix_transform33f3f(I->Matrix,newY,test); printf(" %8.3f %8.3f %8.3f\n",test[0],test[1],test[2]); printf(" %8.3f %8.3f %8.3f\n",I->LightNormal[0],I->LightNormal[1],I->LightNormal[2]); matrix_transform33f3f(I->Matrix,I->LightNormal,test); printf(" %8.3f %8.3f %8.3f\n",test[0],test[1],test[2]); printf(">%8.3f %8.3f %8.3f\n",I->Matrix[0][0],I->Matrix[0][1],I->Matrix[0][2]); printf(">%8.3f %8.3f %8.3f\n",I->Matrix[1][0],I->Matrix[1][1],I->Matrix[1][2]); printf(">%8.3f %8.3f %8.3f\n",I->Matrix[2][0],I->Matrix[2][1],I->Matrix[2][2]); */ } /*========================================================================*/ void BasisGetTriangleFlatDotgle(CBasis *I,RayInfo *r,int i) { float *n0 = I->Normal + (3 * I->Vert2Normal[i]); r->flat_dotgle = n0[2]; } void BasisGetTriangleFlatDotglePerspective(CBasis *I,RayInfo *r,int i) { float *n0 = I->Normal + (3 * I->Vert2Normal[i]); r->flat_dotgle = -dot_product3f(r->dir,n0); } /*========================================================================*/ void BasisGetEllipsoidNormal(CBasis *I,RayInfo *r,int i,int perspective) { if(perspective) { r->impact[0] = r->base[0]+r->dir[0]*r->dist; r->impact[1] = r->base[1]+r->dir[1]*r->dist; r->impact[2] = r->base[2]+r->dir[2]*r->dist; } else { r->impact[0] = r->base[0]; r->impact[1] = r->base[1]; r->impact[2] = r->base[2]-r->dist; } { float *n1 = I->Normal + (3 * I->Vert2Normal[i]); float *n2 = n1 + 3, *n3 = n1 + 6; float *scale = r->prim->n0; float d1,d2,d3,s1,s2,s3; float comp1[3], comp2[3], comp3[3]; float direct[3],surfnormal[3]; direct[0] = r->impact[0]-r->sphere[0]; direct[1] = r->impact[1]-r->sphere[1]; direct[2] = r->impact[2]-r->sphere[2]; normalize3f(direct); d1 = dot_product3f(direct, n1); d2 = dot_product3f(direct, n2); d3 = dot_product3f(direct, n3); if(scale[0]>R_SMALL8) { s1 = d1 / (scale[0] * scale[0]); } else { s1 = 0.0F; } if(scale[1]>R_SMALL8) { s2 = d2 / (scale[1] * scale[1]); } else { s2 = 0.0F; } if(scale[2]>R_SMALL8) { s3 = d3 / (scale[2] * scale[2]); } else { s3 = 0.0F; } scale3f(n1, s1, comp1); scale3f(n2, s2, comp2); scale3f(n3, s3, comp3); copy3f(comp1, surfnormal); add3f(comp2, surfnormal, surfnormal); add3f(comp3, surfnormal, surfnormal); normalize23f(surfnormal, r->surfnormal); } } void BasisGetTriangleNormal(CBasis *I,RayInfo *r,int i,float *fc,int perspective) { float *n0,w2,fc0,fc1,fc2; float vt1[3]; CPrimitive *lprim = r->prim; if(perspective) { r->impact[0] = r->base[0]+r->dir[0]*r->dist; r->impact[1] = r->base[1]+r->dir[1]*r->dist; r->impact[2] = r->base[2]+r->dir[2]*r->dist; } else { r->impact[0] = r->base[0]; r->impact[1] = r->base[1]; r->impact[2] = r->base[2]-r->dist; } n0 = I->Normal + (3 * I->Vert2Normal[i]) + 3; /* skip triangle normal */ w2 = 1.0F - (r->tri1 + r->tri2); /* printf("%8.3f %8.3f\n",r->tri[1],r->tri[2]);*/ fc0 = (lprim->c2[0]*r->tri1)+(lprim->c3[0]*r->tri2)+(lprim->c1[0]*w2); fc1 = (lprim->c2[1]*r->tri1)+(lprim->c3[1]*r->tri2)+(lprim->c1[1]*w2); fc2 = (lprim->c2[2]*r->tri1)+(lprim->c3[2]*r->tri2)+(lprim->c1[2]*w2); r->trans = (lprim->tr[1]*r->tri1)+(lprim->tr[2]*r->tri2)+(lprim->tr[0]*w2); scale3f(n0+3, r->tri1, r->surfnormal); scale3f(n0+6, r->tri2, vt1); add3f(vt1, r->surfnormal, r->surfnormal); scale3f(n0,w2,vt1); add3f(vt1,r->surfnormal,r->surfnormal); normalize3f(r->surfnormal); fc[0] = fc0; fc[1] = fc1; fc[2] = fc2; } /*========================================================================*/ #ifdef PROFILE_BASIS int n_cells = 0; int n_prims = 0; int n_triangles = 0; int n_spheres = 0; int n_cylinders = 0; int n_sausages = 0; int n_skipped = 0; #endif #ifdef _PYMOL_INLINE __inline__ #endif int BasisHitPerspective(BasisCallRec *BC) { register CBasis *BI = BC->Basis; register MapType *map = BI->Map; register int iMin0 = map->iMin[0]; register int iMin1 = map->iMin[1]; register int iMin2 = map->iMin[2]; register int iMax0 = map->iMax[0]; register int iMax1 = map->iMax[1]; register int iMax2 = map->iMax[2]; register int a,b,c; register float iDiv = map->recipDiv; register float base0, base1, base2; register float min0 = map->Min[0] * iDiv; register float min1 = map->Min[1] * iDiv; register float min2 = map->Min[2] * iDiv; int new_ray = !BC->pass; RayInfo *r = BC->rr; MapCache *cache = &BC->cache; int *cache_cache = cache->Cache; int *cache_CacheLink = cache->CacheLink; CPrimitive *r_prim = NULL; if(new_ray) { /* see if we can eliminate this ray right away using the mask */ base0 = (r->base[0] * iDiv) - min0; base1 = (r->base[1] * iDiv) - min1; a = (int)base0; b = (int)base1; a += MapBorder; b += MapBorder; if(a < iMin0) a = iMin0; else if(a > iMax0) a = iMax0; if(b < iMin1) b = iMin1; else if(b > iMax1) b = iMax1; if(!*(map->EMask + a * map->Dim[1] + b)) return -1; } { register int last_a = -1, last_b = -1, last_c = -1; register int allow_break; register int minIndex=-1; register float step0, step1, step2; register float back_dist = BC->back_dist; const float _0 = 0.0F, _1 = 1.0F; float r_tri1=_0, r_tri2=_0, r_dist, dist; /* zero inits to suppress compiler warnings */ float r_sphere0=_0,r_sphere1=_0,r_sphere2=_0; int h,*ip; int excl_trans_flag; int *elist, local_iflag = false; int terminal = -1; int *ehead = map->EHead; int d1d2 = map->D1D2; int d2 = map->Dim[2]; const int *vert2prim = BC->vert2prim; const float excl_trans = BC->excl_trans; const float BasisFudge0 = BC->fudge0; const float BasisFudge1 = BC->fudge1; int v2p; int i,ii; int n_vert = BI->NVertex, n_eElem = map->NEElem; int except1 = BC->except1; int except2 = BC->except2; int check_interior_flag = BC->check_interior; float sph[3],vt[3],tri1=_0,tri2; register CPrimitive *BC_prim = BC->prim; register int *BI_Vert2Normal = BI->Vert2Normal; register float *BI_Vertex = BI->Vertex; register float *BI_Precomp = BI->Precomp; float *BI_Normal = BI->Normal; float *BI_Radius = BI->Radius; float *BI_Radius2 = BI->Radius2; copy3f(r->base, vt); elist = map->EList; r_dist = MAXFLOAT; excl_trans_flag = (excl_trans != _0); if(except1 >= 0) except1 = vert2prim[except1]; if(except2 >= 0) except2 = vert2prim[except2]; MapCacheReset(cache); { /* take steps with a Z-size equil to the grid spacing */ register float div = iDiv * (-MapGetDiv(BI->Map)/r->dir[2]); step0 = r->dir[0]*div; step1 = r->dir[1]*div; step2 = r->dir[2]*div; } base0 = (r->skip[0] * iDiv) - min0; base1 = (r->skip[1] * iDiv) - min1; base2 = (r->skip[2] * iDiv) - min2; allow_break = false; while(1) { int inside_code; register int clamped; a = ((int)base0); b = ((int)base1); c = ((int)base2); inside_code = 1; clamped = false; a += MapBorder; b += MapBorder; c += MapBorder; #define EDGE_ALLOWANCE 1 if(a < iMin0) { if(((iMin0 - a) > EDGE_ALLOWANCE) && allow_break) break; else { a = iMin0; clamped = true; } } else if(a > iMax0) { if(((a - iMax0) > EDGE_ALLOWANCE) && allow_break) break; else { a = iMax0; clamped = true; } } if(b < iMin1) { if(((iMin1 - b) > EDGE_ALLOWANCE) && allow_break) break; else { b = iMin1; clamped = true; } } else if(b > iMax1) { if(((b - iMax1) > EDGE_ALLOWANCE) && allow_break) break; else { b = iMax1; clamped = true; } } if(c < iMin2) { if((iMin2 - c) > EDGE_ALLOWANCE) break; else { c = iMin2; clamped = true; } } else if(c > iMax2) { if((c - iMax2) > EDGE_ALLOWANCE) inside_code = 0; else { c = iMax2; clamped = true; } } if(inside_code && (((a!=last_a)||(b!=last_b)||(c!=last_c)))) { register int new_min_index; h = *(ehead + (a * d1d2) + (b * d2) + c); new_min_index = -1; if(!clamped) /* don't discard a ray until it has hit the objective at least once */ allow_break = true; if((terminal>0)&&(last_c!=c)) { if(!terminal--) break; } if((h>0) && (h<n_eElem)) { int do_loop; ip = elist + h; last_a = a; i = *(ip++); last_b = b; do_loop = ((i>=0)&&(i<n_vert)); last_c = c; while(do_loop) /* n_vert checking is a bug workaround */ { ii = *(ip++); v2p = vert2prim[i]; do_loop = ((ii>=0)&&(ii<n_vert)); if((v2p != except1) && (v2p != except2) && (!MapCached(cache,v2p))) { register CPrimitive *prm = BC_prim + v2p; int prm_type; /*MapCache(cache,v2p);*/ cache_cache[v2p] = 1; prm_type = prm->type; cache_CacheLink[v2p] = cache->CacheStart; cache->CacheStart = v2p; switch(prm_type) { case cPrimTriangle: case cPrimCharacter: { float *dir = r->dir; float *d10 = BI_Precomp + BI_Vert2Normal[i] * 3; float *d20 = d10+3; register float det, inv_det; register float pvec0,pvec1, pvec2; register float dir0 = dir[0], dir1 = dir[1], dir2 = dir[2]; register float d20_0 = d20[0], d20_1 = d20[1], d20_2 = d20[2]; register float d10_0 = d10[0], d10_1 = d10[1], d10_2 = d10[2]; /* cross_product3f(dir, d20, pvec); */ pvec0 = dir1*d20_2 - dir2*d20_1; pvec1 = dir2*d20_0 - dir0*d20_2; pvec2 = dir0*d20_1 - dir1*d20_0; /* det = dot_product3f(pvec, d10); */ det = pvec0 * d10_0 + pvec1 * d10_1 + pvec2 * d10_2; if(fabs(det) >= EPSILON) { float *v0 = BI_Vertex + prm->vert*3; register float tvec0, tvec1, tvec2; register float qvec0, qvec1, qvec2; inv_det = _1/det; /* subtract3f(vt,v0,tvec); */ tvec0 = vt[0] - v0[0]; tvec1 = vt[1] - v0[1]; tvec2 = vt[2] - v0[2]; /* dot_product3f(tvec,pvec) * inv_det; */ tri1 = (tvec0 * pvec0 + tvec1 * pvec1 + tvec2 * pvec2) * inv_det; /* cross_product3f(tvec,d10,qvec); */ qvec0 = tvec1*d10_2 - tvec2*d10_1; qvec1 = tvec2*d10_0 - tvec0*d10_2; if( (tri1>= BasisFudge0) && (tri1<=BasisFudge1)) { qvec2 = tvec0*d10_1 - tvec1*d10_0; /* dot_product3f(dir, qvec) * inv_det; */ tri2 = (dir0 * qvec0 + dir1 * qvec1 + dir2 * qvec2) * inv_det; /* dot_product3f(d20, qvec) * inv_det; */ dist = (d20_0 * qvec0 + d20_1 * qvec1 + d20_2 * qvec2) * inv_det; if( (tri2 >= BasisFudge0) && (tri2 <= BasisFudge1) && ((tri1 + tri2) <= BasisFudge1)) { if( (dist < r_dist) && (dist >= _0) && (dist <= back_dist) && (prm->trans != _1) ) { new_min_index = prm->vert; r_tri1 = tri1; r_tri2 = tri2; r_dist = dist; } } } } } break; case cPrimSphere: { if(LineClipPoint( r->base, r->dir, BI_Vertex + i*3, &dist, BI_Radius[i] , BI_Radius2[i])) { if((dist < r_dist) && (prm->trans != _1)) { if((dist >= _0) && (dist <= back_dist)) { new_min_index = prm->vert; r_dist = dist; } else if(check_interior_flag && (dist<=back_dist)) { if(diffsq3f(vt,BI_Vertex+i*3) < BI_Radius2[i]) { local_iflag = true; r_prim = prm; r_dist = _0; new_min_index = prm->vert; } } } } } break; case cPrimEllipsoid: { if(LineClipPoint( r->base, r->dir, BI_Vertex + i*3, &dist, BI_Radius[i] , BI_Radius2[i])) { if((dist < r_dist) && (prm->trans != _1)) { float *n1 = BI_Normal + BI_Vert2Normal[i] * 3; if(LineClipEllipsoidPoint( r->base, r->dir, BI_Vertex + i*3, &dist, BI_Radius[i] , BI_Radius2[i], prm->n0, n1, n1+3, n1+6)) { if(dist < r_dist) { if((dist >= _0) && (dist <= back_dist)) { new_min_index = prm->vert; r_dist = dist; } else if(check_interior_flag && (dist<=back_dist)) { #if 0 if(diffsq3f(vt,BI_Vertex+i*3) < BI_Radius2[i]) { /* TO FIX */ local_iflag = true; r_prim = prm; r_dist = _0; new_min_index = prm->vert; } #endif } } } } } } break; case cPrimCylinder: if(LineToSphereCapped(r->base,r->dir,BI_Vertex+i*3, BI_Normal+BI_Vert2Normal[i]*3, BI_Radius[i], prm->l1,sph,&tri1, prm->cap1,prm->cap2)) { if(LineClipPoint(r->base,r->dir,sph,&dist,BI_Radius[i],BI_Radius2[i])) { if((dist < r_dist) && (prm->trans != _1)) { if((dist >= _0) && (dist <= back_dist)) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; new_min_index = prm->vert; r_dist = dist; } else if(check_interior_flag && (dist<=back_dist)) { if(FrontToInteriorSphereCapped(vt, BI_Vertex+i*3, BI_Normal+BI_Vert2Normal[i]*3, BI_Radius[i], BI_Radius2[i], prm->l1, prm->cap1, prm->cap2)) { local_iflag = true; r_prim = prm; r_dist = _0; new_min_index = prm->vert; } } } } } break; case cPrimCone: { float sph_rad,sph_rad_sq; if(ConeLineToSphereCapped(r->base,r->dir,BI_Vertex+i*3, BI_Normal+BI_Vert2Normal[i]*3, BI_Radius[i],prm->r2,prm->l1,sph, &tri1, &sph_rad, &sph_rad_sq, prm->cap1,prm->cap2)) { if(LineClipPoint(r->base, r->dir, sph, &dist, sph_rad, sph_rad_sq)) { if((dist < r_dist) && (prm->trans != _1)) { if((dist >= _0) && (dist <= back_dist)) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; /* color blending */ r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; new_min_index = prm->vert; r_dist = dist; } else if(check_interior_flag && (dist<=back_dist)) { if(FrontToInteriorSphereCapped(vt, BI_Vertex+i*3, BI_Normal+BI_Vert2Normal[i]*3, BI_Radius[i], BI_Radius2[i], prm->l1, prm->cap1, prm->cap2)) { local_iflag = true; r_prim = prm; r_dist = _0; new_min_index = prm->vert; } } } } } } break; case cPrimSausage: if(LineToSphere(r->base,r->dir, BI_Vertex+i*3,BI_Normal+BI_Vert2Normal[i]*3, BI_Radius[i],prm->l1,sph,&tri1)) { if(LineClipPoint(r->base,r->dir,sph,&dist,BI_Radius[i],BI_Radius2[i])) { int tmp_flag = false; if((dist<r_dist) && (prm->trans != _1) ) { if((dist >= _0) && (dist <= back_dist)) { tmp_flag = true; if(excl_trans_flag) { if( (prm->trans > _0) && (dist < excl_trans) ) tmp_flag = false; } if(tmp_flag) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; new_min_index = prm->vert; r_dist = dist; } } else if(check_interior_flag &&( dist<=back_dist) ) { if(FrontToInteriorSphere(vt, BI_Vertex+i*3, BI_Normal+BI_Vert2Normal[i]*3, BI_Radius[i], BI_Radius2[i], prm->l1)) { local_iflag = true; r_prim = prm; r_dist = _0; new_min_index = prm->vert; } } } } } break; } /* end of switch */ } /* end of if */ i = ii; } /* end of while */ if(local_iflag) { r->prim = r_prim; r->dist = r_dist; break; } if( new_min_index > -1 ) { minIndex = new_min_index; r_prim = BC_prim + vert2prim[minIndex]; if((r_prim->type == cPrimSphere) || (r_prim->type == cPrimEllipsoid)) { const float *vv = BI->Vertex + minIndex * 3; r_sphere0 = vv[0]; r_sphere1 = vv[1]; r_sphere2 = vv[2]; } BC->interior_flag = local_iflag; r->tri1 = r_tri1; r->tri2 = r_tri2; r->prim = r_prim; r->dist = r_dist; r->sphere[0] = r_sphere0; r->sphere[1] = r_sphere1; r->sphere[2] = r_sphere2; } } /* if -- h valid */ } /* end of if */ if(minIndex > -1) { if(terminal<0) terminal = EDGE_ALLOWANCE+1; } base0+=step0; base1+=step1; base2+=step2; /* advance through the map one block at a time -- note that this is a crappy way to walk through the map... */ } BC->interior_flag = local_iflag; return(minIndex); } } #ifdef _PYMOL_INLINE __inline__ #endif int BasisHitOrthoscopic(BasisCallRec *BC) { const float _0 = 0.0F, _1 = 1.0F; float oppSq,dist=_0,sph[3],vt[3],tri1,tri2; int a,b,c,h,*ip; int excl_trans_flag; int check_interior_flag; int *elist, local_iflag = false; float minusZ[3] = { 0.0F, 0.0F, -1.0F }; CBasis *BI = BC->Basis; RayInfo *r = BC->rr; if( MapInsideXY(BI->Map, r->base, &a, &b, &c) ) { register int minIndex=-1; int v2p; int i,ii; int *xxtmp; int do_loop; int except1 = BC->except1; int except2 = BC->except2; int n_vert = BI->NVertex, n_eElem = BI->Map->NEElem; const int *vert2prim = BC->vert2prim; const float front = BC->front; const float back = BC->back; const float excl_trans = BC->excl_trans; const float BasisFudge0 = BC->fudge0; const float BasisFudge1 = BC->fudge1; MapCache *cache = &BC->cache; float r_tri1=_0, r_tri2=_0, r_dist=_0; /* zero inits to suppress compiler warnings */ float r_sphere0=_0,r_sphere1=_0,r_sphere2=_0; CPrimitive *r_prim = NULL; check_interior_flag = BC->check_interior; /* assumption: always heading in the negative Z direction with our vector... */ vt[0] = r->base[0]; vt[1] = r->base[1]; vt[2] = r->base[2] - front; if(except1 >= 0) except1 = vert2prim[except1]; if(except2 >= 0) except2 = vert2prim[except2]; excl_trans_flag = (excl_trans != _0); r_dist = MAXFLOAT; xxtmp = BI->Map->EHead + (a * BI->Map->D1D2) + (b * BI->Map->Dim[2]) + c; MapCacheReset(cache); elist = BI->Map->EList; while(c >= MapBorder) { h = *xxtmp; if((h>0) && (h<n_eElem)) { ip = elist + h; i = *(ip++); do_loop = ((i>=0)&&(i<n_vert)); while(do_loop) { ii = *(ip++); v2p = vert2prim[i]; do_loop = ((ii>=0)&&(ii<n_vert)); if((v2p != except1) && (v2p != except2) && (!MapCached(cache,v2p))) { CPrimitive *prm = BC->prim + v2p; MapCache(cache,v2p); switch(prm->type) { case cPrimTriangle: case cPrimCharacter: if(!prm->cull) { float *pre = BI->Precomp + BI->Vert2Normal[i] * 3; if( pre[6] ) { float *vert0 = BI->Vertex + prm->vert * 3; float tvec0 = vt[0] - vert0[0]; float tvec1 = vt[1] - vert0[1]; tri1 = (tvec0 * pre[4] - tvec1 * pre[3]) * pre[7]; tri2 = -(tvec0 * pre[1] - tvec1 * pre[0]) * pre[7]; if( !( (tri1 < BasisFudge0) || (tri2 < BasisFudge0) || (tri1 > BasisFudge1) || ((tri1 + tri2) > BasisFudge1) ) ) { dist = (r->base[2] - (tri1*pre[2]) - (tri2*pre[5]) - vert0[2]); if( (dist < r_dist) && (dist >= front) && (dist <= back) && (prm->trans != _1) ) { minIndex = prm->vert; r_tri1 = tri1; r_tri2 = tri2; r_dist = dist; } } } } break; case cPrimSphere: oppSq = ZLineClipPoint( r->base, BI->Vertex + i*3, &dist, BI->Radius[i] ); if(oppSq <= BI->Radius2[i]) { dist = (float)(sqrt1f(dist) - sqrt1f((BI->Radius2[i]-oppSq))); if((dist < r_dist) && (prm->trans != _1)) { if((dist >= front) && (dist <= back)) { minIndex = prm->vert; r_dist = dist; } else if(check_interior_flag) { if(diffsq3f(vt,BI->Vertex+i*3) < BI->Radius2[i]) { local_iflag = true; r_prim = prm; r_dist = front; minIndex = prm->vert; } } } } break; case cPrimEllipsoid: oppSq = ZLineClipPoint( r->base, BI->Vertex + i*3, &dist, BI->Radius[i] ); if(oppSq <= BI->Radius2[i]) { dist = (float)(sqrt1f(dist) - sqrt1f((BI->Radius2[i]-oppSq))); if((dist < r_dist) && (prm->trans != _1)) { float *n1 = BI->Normal + BI->Vert2Normal[i] * 3; if(LineClipEllipsoidPoint( r->base, minusZ, BI->Vertex + i*3, &dist, BI->Radius[i] , BI->Radius2[i], prm->n0, n1, n1+3, n1+6)) { if(dist < r_dist) { if((dist >= _0) && (dist <= back)) { minIndex = prm->vert; r_dist = dist; } else if(check_interior_flag) { #if 0 /* TODO TO FIX */ if(diffsq3f(vt,BI->Vertex+i*3) < BI->Radius2[i]) { local_iflag = true; r_prim = prm; r_dist = front; minIndex = prm->vert; } #endif } } } } } break; case cPrimCylinder: if(ZLineToSphereCapped(r->base,BI->Vertex+i*3, BI->Normal+BI->Vert2Normal[i]*3, BI->Radius[i], prm->l1,sph,&tri1,prm->cap1,prm->cap2, BI->Precomp + BI->Vert2Normal[i] * 3)) { oppSq = ZLineClipPoint(r->base,sph,&dist,BI->Radius[i]); if(oppSq<=BI->Radius2[i]) { dist = (float)(sqrt1f(dist)-sqrt1f((BI->Radius2[i]-oppSq))); if((dist < r_dist) && (prm->trans != _1)) { if((dist >= front) && (dist <= back)) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; minIndex = prm->vert; r_dist = dist; } else if(check_interior_flag) { if(FrontToInteriorSphereCapped(vt, BI->Vertex+i*3, BI->Normal+BI->Vert2Normal[i]*3, BI->Radius[i], BI->Radius2[i], prm->l1, prm->cap1, prm->cap2)) { local_iflag = true; r_prim = prm; r_dist = front; minIndex = prm->vert; } } } } } break; case cPrimCone: { float sph_rad,sph_rad_sq; if(ConeLineToSphereCapped(r->base,minusZ,BI->Vertex+i*3, BI->Normal+BI->Vert2Normal[i]*3, BI->Radius[i],prm->r2,prm->l1,sph, &tri1, &sph_rad, &sph_rad_sq, prm->cap1,prm->cap2)) { oppSq = ZLineClipPoint(r->base,sph,&dist,sph_rad); if(oppSq<=sph_rad_sq) { dist = (float)(sqrt1f(dist)-sqrt1f((sph_rad_sq-oppSq))); if((dist < r_dist) && (prm->trans != _1)) { if((dist >= front) && (dist <= back)) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; minIndex = prm->vert; r_dist = dist; } else if(check_interior_flag) { if(FrontToInteriorSphereCapped(vt, BI->Vertex+i*3, BI->Normal+BI->Vert2Normal[i]*3, sph_rad, sph_rad_sq, prm->l1, prm->cap1, prm->cap2)) { local_iflag = true; r_prim = prm; r_dist = front; minIndex = prm->vert; } } } } } } break; case cPrimSausage: if(ZLineToSphere(r->base,BI->Vertex+i*3,BI->Normal+BI->Vert2Normal[i]*3, BI->Radius[i],prm->l1,sph,&tri1, BI->Precomp + BI->Vert2Normal[i] * 3)) { oppSq = ZLineClipPoint(r->base,sph,&dist,BI->Radius[i]); if(oppSq <= BI->Radius2[i]) { int tmp_flag = false; dist = (float)(sqrt1f(dist)-sqrt1f((BI->Radius2[i]-oppSq))); if((dist<r_dist) && (prm->trans != _1) ) { if((dist >= front) && (dist <= back)) { tmp_flag = true; if(excl_trans_flag) { if( (prm->trans > _0) && (dist < excl_trans) ) tmp_flag = false; } if(tmp_flag) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; minIndex = prm->vert; r_dist = dist; } } else if(check_interior_flag) { if(FrontToInteriorSphere(vt, BI->Vertex+i*3, BI->Normal+BI->Vert2Normal[i]*3, BI->Radius[i], BI->Radius2[i], prm->l1)) { local_iflag = true; r_prim = prm; r_dist = front; minIndex = prm->vert; } } } } } break; } /* end of switch */ } /* end of if */ i = ii; } /* end of while */ } /* and of course stop when we hit the edge of the map */ if(local_iflag) break; /* we've processed all primitives associated with this voxel, so if an intersection has been found which occurs in front of the next voxel, then we can stop */ if( minIndex > -1 ) { int aa,bb,cc; vt[2] = r->base[2] - r_dist; MapLocus(BI->Map,vt,&aa,&bb,&cc); if(cc > c) break; else vt[2] = r->base[2] - front; } c--; xxtmp--; } /* end of while */ if( minIndex > -1 ) { r_prim = BC->prim + vert2prim[minIndex]; if((r_prim->type == cPrimSphere) || (r_prim->type == cPrimEllipsoid)) { const float *vv = BI->Vertex + minIndex * 3; r_sphere0 = vv[0]; r_sphere1 = vv[1]; r_sphere2 = vv[2]; } } BC->interior_flag = local_iflag; r->tri1 = r_tri1; r->tri2 = r_tri2; r->prim = r_prim; r->dist = r_dist; r->sphere[0] = r_sphere0; r->sphere[1] = r_sphere1; r->sphere[2] = r_sphere2; return(minIndex); } /* end of if */ BC->interior_flag = local_iflag; return(-1); } #ifdef _PYMOL_INLINE __inline__ #endif int BasisHitShadow(BasisCallRec *BC) { const float _0 = 0.0F; const float _1 = 1.0F; float oppSq,dist=_0,tri1,tri2; float sph[3],vt[3]; register int h,*ip; int a,b,c; int excl_trans_flag; int check_interior_flag; int *elist, local_iflag = false; float minusZ[3] = { 0.0F, 0.0F, -1.0F }; /* local copies (eliminate these extra copies later on) */ CBasis *BI = BC->Basis; RayInfo *r = BC->rr; if( MapInsideXY(BI->Map,r->base, &a, &b, &c) ) { register int minIndex=-1; register int v2p; register int i,ii; int *xxtmp; int n_vert = BI->NVertex, n_eElem = BI->Map->NEElem; register int except1 = BC->except1; register int except2 = BC->except2; const int *vert2prim = BC->vert2prim; const int trans_shadows = BC->trans_shadows; const int nearest_shadow = BC->nearest_shadow; const float excl_trans = BC->excl_trans; const float BasisFudge0 = BC->fudge0; const float BasisFudge1 = BC->fudge1; const int label_shadow_mode = BC->label_shadow_mode; register MapCache *cache = &BC->cache; int *cache_cache = cache->Cache; int *cache_CacheLink = cache->CacheLink; register CPrimitive *BC_prim = BC->prim; register float r_tri1=_0, r_tri2=_0, r_dist; /* zero inits to suppress compiler warnings */ register float r_sphere0=_0,r_sphere1=_0,r_sphere2=_0; register float r_trans = _0; CPrimitive *r_prim = NULL; check_interior_flag = BC->check_interior; /* assumption: always heading in the negative Z direction with our vector... */ vt[0] = r->base[0]; vt[1] = r->base[1]; if(except1 >= 0) except1 = vert2prim[except1]; if(except2 >= 0) except2 = vert2prim[except2]; excl_trans_flag = (excl_trans != _0); r_trans = _1; r_dist = MAXFLOAT; xxtmp = BI->Map->EHead + (a * BI->Map->D1D2) + (b * BI->Map->Dim[2]) + c; MapCacheReset(cache); elist = BI->Map->EList; while(c >= MapBorder) { h = *xxtmp; if((h>0) && (h<n_eElem)) { int do_loop; ip = elist + h; i = *(ip++); do_loop = ((i>=0)&&(i<n_vert)); while(do_loop) { ii = *(ip++); v2p = vert2prim[i]; do_loop = ((ii>=0)&&(ii<n_vert)); if( (v2p != except1) && (v2p != except2) && !MapCached(cache,v2p) ) { register CPrimitive *prm = BC_prim + v2p; int prm_type; /*MapCache(cache,v2p);*/ cache_cache[v2p] = 1; prm_type = prm->type; cache_CacheLink[v2p] = cache->CacheStart; cache->CacheStart = v2p; switch(prm_type) { case cPrimCharacter: /* will need special handling for character shadows */ if(label_shadow_mode&0x2) { /* if labels case shadows... */ float *pre = BI->Precomp + BI->Vert2Normal[i] * 3; if( pre[6] ) { float *vert0 = BI->Vertex + prm->vert * 3; float tvec0 = vt[0] - vert0[0]; float tvec1 = vt[1] - vert0[1]; tri1 = (tvec0 * pre[4] - tvec1 * pre[3]) * pre[7]; tri2 = -(tvec0 * pre[1] - tvec1 * pre[0]) * pre[7]; if( !( (tri1 < BasisFudge0) || (tri2 < BasisFudge0) || (tri1 > BasisFudge1) || ((tri1 + tri2) > BasisFudge1) ) ) { dist = (r->base[2] - (tri1*pre[2]) - (tri2*pre[5]) - vert0[2]); { float fc[3]; float trans; r->tri1 = tri1; r->tri2 = tri2; r->dist = dist; r->prim = prm; { float w2; w2 = _1 - (r->tri1 + r->tri2); fc[0] = (prm->c2[0]*r->tri1)+(prm->c3[0]*r->tri2)+(prm->c1[0]*w2); fc[1] = (prm->c2[1]*r->tri1)+(prm->c3[1]*r->tri2)+(prm->c1[1]*w2); fc[2] = (prm->c2[2]*r->tri1)+(prm->c3[2]*r->tri2)+(prm->c1[2]*w2); } trans = CharacterInterpolate(BI->G,prm->char_id,fc); if(trans == _0) { /* opaque? return immed. */ if(dist > -kR_SMALL4) { if(nearest_shadow) { if(dist < r_dist) { minIndex = prm->vert; r_tri1 = tri1; r_tri2 = tri2; r_dist = dist; r_trans = (r->trans = trans); } } else { r->prim = prm; r->trans = _0; r->dist = dist; return(1); } } } else if(trans_shadows) { if((dist > -kR_SMALL4) && (( r_trans > trans) || (nearest_shadow && (dist<r_dist) && ( r_trans >= trans)))) { minIndex = prm->vert; r_tri1 = tri1; r_tri2 = tri2; r_dist = dist; r_trans = (r->trans = trans); } } } } } } break; case cPrimTriangle: { float *pre = BI->Precomp + BI->Vert2Normal[i] * 3; if( pre[6] ) { float *vert0 = BI->Vertex + prm->vert * 3; float tvec0 = vt[0] - vert0[0]; float tvec1 = vt[1] - vert0[1]; tri1 = (tvec0 * pre[4] - tvec1 * pre[3]) * pre[7]; tri2 = -(tvec0 * pre[1] - tvec1 * pre[0]) * pre[7]; if( !( (tri1 < BasisFudge0) || (tri2 < BasisFudge0) || (tri1 > BasisFudge1) || ((tri1+tri2)> BasisFudge1) ) ) { float *tr = prm->tr; float trans = _0; dist = (r->base[2] - (tri1*pre[2]) - (tri2*pre[5]) - vert0[2]); if(prm->trans!=_0) { trans = (tr[1]*tri1)+(tr[2]*tri2)+(tr[0]*(_1 - (tri1+tri2))); } if(trans == _0) { if(dist > -kR_SMALL4) { if(nearest_shadow) { /* do we need the nearest shadow? */ if(dist < r_dist) { minIndex = prm->vert; r_tri1 = tri1; r_tri2 = tri2; r_dist = dist; r_trans = (r->trans = trans); } } else { r->prim = prm; r->trans = _0; r->dist = dist; return(1); } } } else if(trans_shadows) { if((dist > -kR_SMALL4) && (( r_trans > trans) || (nearest_shadow && (dist<r_dist) && ( r_trans >= trans)))) { minIndex = prm->vert; r_tri1 = tri1; r_tri2 = tri2; r_dist = dist; r_trans = (r->trans = trans); } } } } } break; case cPrimSphere: oppSq = ZLineClipPoint( r->base, BI->Vertex + i*3, &dist, BI->Radius[i] ); if(oppSq <= BI->Radius2[i]) { dist = (float)(sqrt1f(dist) - sqrt1f((BI->Radius2[i]-oppSq))); if(prm->trans == _0) { if(dist > -kR_SMALL4) { if(nearest_shadow) { if(dist < r_dist) { minIndex = prm->vert; r_dist = dist; r_trans = (r->trans = prm->trans); } } else { r->prim = prm; r->trans = prm->trans; r->dist = dist; return(1); } } } else if(trans_shadows) { if((dist > -kR_SMALL4) && (( r_trans > prm->trans) || (nearest_shadow && (dist<r_dist) && ( r_trans >= prm->trans)))) { minIndex = prm->vert; r_dist = dist; r_trans = (r->trans = prm->trans); } } } break; case cPrimEllipsoid: oppSq = ZLineClipPointNoZCheck( r->base, BI->Vertex + i*3, &dist, BI->Radius[i] ); if(oppSq <= BI->Radius2[i]) { dist = (float)(sqrt1f(dist) - sqrt1f((BI->Radius2[i]-oppSq))); if((dist < r_dist ) || (trans_shadows && (r_trans!=_0))) { float *n1 = BI->Normal + BI->Vert2Normal[i] * 3; if(LineClipEllipsoidPoint( r->base, minusZ, BI->Vertex + i*3, &dist, BI->Radius[i] , BI->Radius2[i], prm->n0, n1, n1+3, n1+6)) { if(prm->trans == _0) { if(dist > -kR_SMALL4) { if(nearest_shadow) { if(dist < r_dist) { minIndex = prm->vert; r_dist = dist; r_trans = (r->trans = prm->trans); } } else { r->prim = prm; r->trans = prm->trans; r->dist = dist; return(1); } } } else if(trans_shadows) { if((dist > -kR_SMALL4) && (( r_trans > prm->trans) || (nearest_shadow && (dist<r_dist) && ( r_trans >= prm->trans)))) { minIndex = prm->vert; r_dist = dist; r_trans = (r->trans = prm->trans); } } } } } break; case cPrimCone: { float sph_rad,sph_rad_sq; if(ConeLineToSphereCapped(r->base,minusZ,BI->Vertex+i*3, BI->Normal+BI->Vert2Normal[i]*3, BI->Radius[i],prm->r2,prm->l1,sph, &tri1, &sph_rad, &sph_rad_sq, 1,1)) { oppSq = ZLineClipPoint(r->base,sph,&dist,sph_rad); if(oppSq <= sph_rad_sq) { dist=(float)(sqrt1f(dist)-sqrt1f((sph_rad_sq-oppSq))); if(prm->trans == _0) { if(dist > -kR_SMALL4) { if(nearest_shadow) { if(dist < r_dist) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; minIndex = prm->vert; r->trans = prm->trans; r_dist = dist; r_trans = (r->trans = prm->trans);} } else { r->prim = prm; r->trans = prm->trans; r->dist = dist; return(1); } } } else if(trans_shadows) { if((dist > -kR_SMALL4) && (( r_trans > prm->trans) || (nearest_shadow && (dist<r_dist) && ( r_trans >= prm->trans)))) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; minIndex = prm->vert; r->trans = prm->trans; r_dist = dist; r_trans = (r->trans = prm->trans); } } } } } break; case cPrimCylinder: if(ZLineToSphereCapped(r->base,BI->Vertex+i*3, BI->Normal+BI->Vert2Normal[i]*3, BI->Radius[i], prm->l1,sph,&tri1,prm->cap1,prm->cap2, BI->Precomp + BI->Vert2Normal[i] * 3)) { oppSq = ZLineClipPoint(r->base,sph,&dist,BI->Radius[i]); if(oppSq <= BI->Radius2[i]) { dist=(float)(sqrt1f(dist)-sqrt1f((BI->Radius2[i]-oppSq))); if(prm->trans == _0) { if(dist > -kR_SMALL4) { if(nearest_shadow) { if(dist < r_dist) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; minIndex = prm->vert; r->trans = prm->trans; r_dist = dist; r_trans = (r->trans = prm->trans);} } else { r->prim = prm; r->trans = prm->trans; r->dist = dist; return(1); } } } else if(trans_shadows) { if((dist > -kR_SMALL4) && (( r_trans > prm->trans) || (nearest_shadow && (dist<r_dist) && ( r_trans >= prm->trans)))) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; minIndex = prm->vert; r->trans = prm->trans; r_dist = dist; r_trans = (r->trans = prm->trans); } } } } break; case cPrimSausage: if(ZLineToSphere(r->base,BI->Vertex+i*3,BI->Normal+BI->Vert2Normal[i]*3, BI->Radius[i],prm->l1,sph,&tri1, BI->Precomp + BI->Vert2Normal[i] * 3)) { oppSq = ZLineClipPoint(r->base,sph,&dist,BI->Radius[i]); if(oppSq <= BI->Radius2[i]) { dist = (float)(sqrt1f(dist) - sqrt1f((BI->Radius2[i]-oppSq))); if(prm->trans == _0) { if(dist > -kR_SMALL4) { if(nearest_shadow) { if(dist < r_dist) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; minIndex = prm->vert; r_dist = dist; r_trans = (r->trans = prm->trans); } } else { r->prim = prm; r->trans = prm->trans; r->dist = dist; return(1); } } } else if(trans_shadows) { if((dist > -kR_SMALL4) && (( r_trans > prm->trans) || (nearest_shadow && (dist<r_dist) && ( r_trans >= prm->trans)))) { if(prm->l1 > kR_SMALL4) r_tri1 = tri1 / prm->l1; r_sphere0 = sph[0]; r_sphere1 = sph[1]; r_sphere2 = sph[2]; minIndex = prm->vert; r_dist = dist; r_trans = (r->trans = prm->trans); } } } } break; } /* end of switch */ } /* end of if */ i = ii; } /* end of while */ } /* and of course stop when we hit the edge of the map */ if(local_iflag) break; /* we've processed all primitives associated with this voxel, so if an intersection has been found which occurs in front of the next voxel, then we can stop */ /* this optimization invalid for transparent surfaces if( minIndex > -1 ) { int aa,bb,cc; vt[2] = r->base[2] - r_dist; MapLocus(BI->Map,vt,&aa,&bb,&cc); if(cc > c) break; } */ c--; xxtmp--; } /* end of while */ if( minIndex > -1 ) { r_prim = BC->prim + vert2prim[minIndex]; if((r_prim->type == cPrimSphere) || (r_prim->type == cPrimEllipsoid)) { const float *vv = BI->Vertex + minIndex * 3; r_sphere0 = vv[0]; r_sphere1 = vv[1]; r_sphere2 = vv[2]; } } BC->interior_flag = local_iflag; r->tri1 = r_tri1; r->tri2 = r_tri2; r->prim = r_prim; r->dist = r_dist; r->trans = r_trans; r->sphere[0] = r_sphere0; r->sphere[1] = r_sphere1; r->sphere[2] = r_sphere2; return(minIndex); } /* end of if */ BC->interior_flag = local_iflag; return(-1); } #if 0 /*========================================================================*/ void BasisOptimizeMap(CBasis *I,float *vertex,int n,int *vert2prim); void BasisOptimizeMap(CBasis *I,float *vertex,int n,int *vert2prim) { int *ip,*op; int k; int v2p; int i,ii; MapCache cache; int a,b,c,aa,bb,cc; int h; int *xxtmp,*elist; MapCacheInit(&cache,I->Map); int cnt = 0,tcnt = 0; elist = I->Map->EList; #if 0 for(k=0;k<n;k++) { if (MapExclLocus(I->Map,vertex + 3*k,&aa,&bb,&cc)) { for (a=aa-1;a<=aa+1;a++) { for (b=bb-1;b<=bb+1;b++) { xxtmp = I->Map->EHead + (a * I->Map->D1D2) + (b * I->Map->Dim[2]); for(c=cc-1;c<=cc+1;c++) { #else { { for( a = I->Map->iMin[0]; a<=I->Map->iMax[0]; a++) { for(b = I->Map->iMin[1]; b<=I->Map->iMax[1]; b++) { xxtmp = I->Map->EHead + (a * I->Map->D1D2) + (b * I->Map->Dim[2]); for( c = I->Map->iMin[2]; c <= I->Map->iMax[2]; c++) { #endif h = *(xxtmp + c); cnt = 0; tcnt = 0; if(h) { op = ip = elist + h; MapCacheReset(&cache); if(ip) { i = *(ip++); while(i >= 0) { tcnt++; v2p = vert2prim[i]; ii = *(ip++); /* printf("%d\n",v2p);*/ if( !MapCached(&cache,v2p)) { *(op++) = i; /* copy/store */ MapCache(&cache,v2p); } /* end of if */ else { /* remove this vertex from the linked list -- we don't need it more than once */ cnt++; } i = ii; } /* end of while */ *op = -1; /* terminate new list */ } /* printf("%6d %6d\n",tcnt,cnt);*/ } } } } } } } #endif /*========================================================================*/ void BasisMakeMap(CBasis *I,int *vert2prim,CPrimitive *prim,int n_prim, float *volume, int group_id,int block_base, int perspective,float front,float size_hint) { register float *v; float ll; CPrimitive *prm; register int i; register int *tempRef; int n,h,q,x,y,z,j,k,l,e; int extra_vert = 0; float p[3],dd[3],*d1,*d2,vd[3],cx[3],cy[3]; float *tempVertex; float xs,ys; int remapMode = true; /* remap mode means that some objects will span more * than one voxel, so we have to worry about populating * those voxels and also about eliminating duplicates * when traversing the neighbor lists */ float min[3],max[3],extent[6]; float sep; float diagonal[3]; float l1,l2; float bh,ch; int n_voxel; const float _0 = 0.0; const float _p5 = 0.5; PRINTFD(I->G,FB_Ray) " BasisMakeMap: I->NVertex %d [(%8.3f, %8.3f, %8.3f),...]\n",I->NVertex, I->Vertex[0],I->Vertex[1],I->Vertex[2] ENDFD; sep = I->MinVoxel; if(sep == _0) { remapMode = false; sep = I->MaxRadius; /* also will imply no remapping of vertices */ } /* we need to get a sense of the actual size in order to avoid sep being too small */ v = I->Vertex; min[0] = max[0] = v[0]; min[1] = max[1] = v[1]; min[2] = max[2] = v[2]; v += 3; { register int a; for(a = 1; a < I->NVertex; a++) { if(min[0] > v[0]) min[0] = v[0]; if(max[0] < v[0]) max[0] = v[0]; if(min[1] > v[1]) min[1] = v[1]; if(max[1] < v[1]) max[1] = v[1]; if(min[2] > v[2]) min[2] = v[2]; if(max[2] < v[2]) max[2] = v[2]; v += 3; } } if(volume) { /* if(min[0] > volume[0]) min[0] = volume[0]; if(max[0] < volume[1]) max[0] = volume[1]; if(min[1] > volume[2]) min[1] = volume[2]; if(max[1] < volume[3]) max[1] = volume[3]; if(min[2] > (-volume[5])) min[2] = (-volume[5]); if(max[2] < (-volume[4])) max[2] = (-volume[4]); */ if(min[0] < volume[0]) min[0] = volume[0]; if(max[0] > volume[1]) max[0] = volume[1]; if(min[1] < volume[2]) min[1] = volume[2]; if(max[1] > volume[3]) max[1] = volume[3]; if(min[2] > (-volume[5])) min[2] = (-volume[5]); if(max[2] < (-volume[4])) max[2] = (-volume[4]); PRINTFB(I->G,FB_Ray,FB_Debugging) " BasisMakeMap: (%8.3f,%8.3f),(%8.3f,%8.3f),(%8.3f,%8.3f)\n", volume[0],volume[1],volume[2],volume[3],volume[4],volume[5] ENDFB(I->G); } /* don't break up space unnecessarily if we only have a few vertices... */ if(I->NVertex) { l1 = (float)fabs(max[0]-min[0]); l2 = (float)fabs(max[1]-min[1]); if(l1 < l2) l1 = l2; l2 = (float)fabs(max[2]-min[2]); if(l1 < l2) l1 = l2; if(l1 < kR_SMALL4) l1 = 100.0; if(I->NVertex < (l1/sep)) sep = (l1/I->NVertex); } if(Feedback(I->G,FB_Ray,FB_Debugging)) { dump3f(min," BasisMakeMap: min"); dump3f(max," BasisMakeMap: max"); dump3f(I->Vertex," BasisMakeMap: I->Vertex"); fflush(stdout); } sep = MapGetSeparation(I->G,sep,max,min,diagonal); /* this needs to be a minimum * estimate of the actual value */ /* here we have to carry out a complicated work-around in order to * efficiently encode our primitives into the map in a way that doesn't * require expanding the map cutoff to the size of the largest object*/ if(remapMode) { register int a,b,c; if(sep<size_hint) /* this keeps us from wasting time & memory on unnecessary subdivision */ sep = size_hint; for(a = 0; a < I->NVertex; a++) { prm = prim + vert2prim[a]; switch(prm->type) { case cPrimTriangle: case cPrimCharacter: if(a == prm->vert) { /* only do this calculation for one of the three vertices */ l1 = (float)length3f(I->Precomp+I->Vert2Normal[a]*3); l2 = (float)length3f(I->Precomp+I->Vert2Normal[a]*3+3); if((l1>=sep)||(l2>=sep)) { b = (int)ceil(l1/sep)+1; c = (int)ceil(l2/sep)+1; extra_vert += 4*b*c; } } break; case cPrimCone: case cPrimCylinder: case cPrimSausage: if((prm->l1+2*prm->r1)>=sep) { q = ((int)(2*(floor(prm->r1/sep)+1)))+1; q = q * q * ((int)ceil((prm->l1+2*prm->r1)/sep)+2); extra_vert+= q; } break; case cPrimEllipsoid: case cPrimSphere: if(prm->r1>=sep) { b = (int)(2*floor(prm->r1/sep)+1); extra_vert += (b*b*b); } break; } } /* for */ extra_vert += I->NVertex; tempVertex = CacheAlloc(I->G,float,extra_vert*3,group_id,cCache_basis_tempVertex); tempRef = CacheAlloc(I->G,int,extra_vert,group_id,cCache_basis_tempRef); ErrChkPtr(I->G,tempVertex); /* can happen if extra vert is unreasonable */ ErrChkPtr(I->G,tempRef); /* lower indexes->flags, top is ref->lower index*/ { register float *vv,*d; n = I->NVertex; v = tempVertex; vv = I->Vertex; memcpy( v, vv, n * sizeof(float) * 3 ); vv += n * 3; v += n * 3; for(a = 0; a < I->NVertex; a++) { prm = prim + vert2prim[a]; switch(prm->type) { case cPrimTriangle: case cPrimCharacter: if(a == prm->vert) { /* only do this calculation for one of the three vertices */ d1 = I->Precomp + I->Vert2Normal[a]*3; d2 = I->Precomp + I->Vert2Normal[a]*3+3; vv = I->Vertex + a*3; l1 = (float)length3f(d1); l2 = (float)length3f(d2); if((l1>=sep)||(l2>=sep)) { b = (int)floor(l1/sep)+1; c = (int)floor(l2/sep)+1; extra_vert += b*c; bh = (float)(b/2)+1; ch = (float)(c/2)+1; for(x = 0; x < bh; x++) { const float xb = (float)x / b; for(y = 0; y < ch; y++) { const float yc = (float)y / c; *(v++) = vv[0] + (d1[0] * xb) + (d2[0] * yc); *(v++) = vv[1] + (d1[1] * xb) + (d2[1] * yc); *(v++) = vv[2] + (d1[2] * xb) + (d2[2] * yc); tempRef[n++] = a; } } for(x = 0; x < bh; x++) { const float xb = (float)x / b; for(y = 0; y < ch; y++) { const float yc = (float)y/c; if( (xb + yc) < _p5 ) { *(v++) = vv[0] + d1[0] * (_p5 + xb) + (d2[0]*yc); *(v++) = vv[1] + d1[1] * (_p5 + xb) + (d2[1]*yc); *(v++) = vv[2] + d1[2] * (_p5 + xb) + (d2[2]*yc); tempRef[n++] = a; *(v++) = vv[0] + (d1[0]*xb) + d2[0]*(_p5 + yc); *(v++) = vv[1] + (d1[1]*xb) + d2[1]*(_p5 + yc); *(v++) = vv[2] + (d1[2]*xb) + d2[2]*(_p5 + yc); tempRef[n++] = a; } } } } } /* if */ break; case cPrimCone: case cPrimCylinder: case cPrimSausage: ll = prm->l1+2*prm->r1; if(ll>=sep) { d = I->Normal+3*I->Vert2Normal[a]; vv = I->Vertex+a*3; get_system1f3f(d,cx,cy); /* creates an orthogonal system about d */ p[0] = vv[0]-d[0]*prm->r1; p[1] = vv[1]-d[1]*prm->r1; p[2] = vv[2]-d[2]*prm->r1; dd[0] = d[0]*sep; dd[1] = d[1]*sep; dd[2] = d[2]*sep; q = (int)floor(prm->r1 / sep) + 1; while(1) { vd[0] = (p[0] += dd[0]); vd[1] = (p[1] += dd[1]); vd[2] = (p[2] += dd[2]); for(x = -q; x <= q; x++) { for(y = -q; y <= q; y++) { xs = x*sep; ys = y*sep; *(v++) = vd[0] + xs*cx[0] + ys*cy[0]; *(v++) = vd[1] + xs*cx[1] + ys*cy[1]; *(v++) = vd[2] + xs*cx[2] + ys*cy[2]; tempRef[n++] = a; } } if(ll <= _0) break; ll -= sep; } } break; case cPrimEllipsoid: case cPrimSphere: if(prm->r1>=sep) { q = (int)floor(prm->r1 / sep); vv = I->Vertex + a*3; for(x = -q; x <= q; x++) { for(y = -q; y <= q; y++) { for( z = -q; z <= q; z++) { *(v++) = vv[0] + x * sep; *(v++) = vv[1] + y * sep; *(v++) = vv[2] + z * sep; tempRef[n++] = a; } } } } break; } /* end of switch */ } } if(n > extra_vert) { printf("BasisMakeMap: %d>%d\n",n,extra_vert); ErrFatal(I->G,"BasisMakeMap","used too many extra vertices (this indicates a bug)...\n"); } PRINTFB(I->G,FB_Ray,FB_Blather) " BasisMakeMap: %d total vertices\n",n ENDFB(I->G); if(volume) { v = tempVertex; min[0] = max[0] = v[0]; min[1] = max[1] = v[1]; min[2] = max[2] = v[2]; v += 3; if(Feedback(I->G,FB_Ray,FB_Debugging)) { dump3f(min," BasisMakeMap: remapped min"); dump3f(max," BasisMakeMap: remapped max"); fflush(stdout); } for(a = 1; a < n; a++) { if(min[0] > v[0]) min[0] = v[0]; if(max[0] < v[0]) max[0] = v[0]; if(min[1] > v[1]) min[1] = v[1]; if(max[1] < v[1]) max[1] = v[1]; if(min[2] > v[2]) min[2] = v[2]; if(max[2] < v[2]) max[2] = v[2]; v += 3; } if(Feedback(I->G,FB_Ray,FB_Debugging)) { dump3f(min," BasisMakeMap: remapped min"); dump3f(max," BasisMakeMap: remapped max"); fflush(stdout); } if(min[0] < volume[0]) min[0] = volume[0]; if(max[0] > volume[1]) max[0] = volume[1]; if(min[1] < volume[2]) min[1] = volume[2]; if(max[1] > volume[3]) max[1] = volume[3]; if(min[2] < (-volume[5])) min[2] = (-volume[5]); if(max[2] > (-volume[4])) max[2] = (-volume[4]); extent[0] = min[0]; extent[1] = max[0]; extent[2] = min[1]; extent[3] = max[1]; extent[4] = min[2]; extent[5] = max[2]; PRINTFB(I->G,FB_Ray,FB_Blather) " BasisMakeMap: Extent [%8.2f %8.2f] [%8.2f %8.2f] [%8.2f %8.2f]\n", extent[0],extent[1],extent[2],extent[3],extent[4],extent[5] ENDFB(I->G); I->Map = MapNewCached(I->G,-sep,tempVertex,n,extent,group_id,block_base); } else { I->Map = MapNewCached(I->G,sep,tempVertex,n,NULL,group_id,block_base); } n_voxel = I->Map->Dim[0]*I->Map->Dim[1]*I->Map->Dim[2]; if(perspective) { /* this is a new optimization which prevents primitives contained entirely within a single voxel from spilling over into neighboring voxels for performance of intersection checks (NOTE: this currently only helps with triangles, and is only useful in optimizing between Z layers). The main impact of this optimization is to reduce the size of the express list (I->Map->Elist) array by about 3-fold */ int *prm_spanner = Calloc(int,n_prim); int *spanner = Calloc(int,n); float *v=tempVertex; int j,k,l,jj,kk,ll; int nVertex = I->NVertex; int prm_index; MapType *map = I->Map; /* figure out which primitives span more than one voxel */ for(a=0;a<n;a++) { if(a<nVertex) prm_index = vert2prim[a]; else prm_index = vert2prim[tempRef[a]]; if( !prm_spanner[prm_index] ) { prm = prim + prm_index; { float *vv = prm->vert*3 + I->Vertex; switch(prm->type) { case cPrimTriangle: case cPrimCharacter: MapLocus(map, v, &j, &k, &l ); MapLocus(map, vv, &jj, &kk, &ll); if((j!=jj)||(k!=kk)||(l!=ll)) { prm_spanner[prm_index] = 1; #if 0 zero3f(prm->c1); zero3f(prm->c2); zero3f(prm->c3); prm->c1[0]=1.0F; prm->c2[0]=1.0F; prm->c3[0]=1.0F; #endif } break; default: /* currently we aren't optimizing other primitives */ prm_spanner[prm_index] = 1; break; } } } v+=3; } /* now flag associated vertices */ for(a=0;a<n;a++) { if(a<nVertex) prm_index = vert2prim[a]; else prm_index = vert2prim[tempRef[a]]; spanner[a] = prm_spanner[prm_index]; } /* and do the optimized expansion */ MapSetupExpressPerp(I->Map,tempVertex,front,n,true,spanner); FreeP(spanner); FreeP(prm_spanner); } else if(n_voxel < (3*n)) { MapSetupExpressXY(I->Map,n,true); } else { MapSetupExpressXYVert(I->Map,tempVertex,n,true); } { register MapType *map = I->Map; register int *sp,*ip,*ip0,ii; register int *elist = map->EList, *ehead = map->EHead; register int *elist_new = elist, *ehead_new = ehead; register int newelem = 0, neelem = -map->NEElem; register int i_nVertex = I->NVertex; const int iMin0=map->iMin[0]; const int iMin1=map->iMin[1]; const int iMin2=map->iMin[2]; const int iMax0=map->iMax[0]; const int iMax1=map->iMax[1]; const int iMax2=map->iMax[2]; #if 0 elist_new = VLACacheAlloc(I->G,int,map->NEElem,group_id,block_base + cCache_map_elist_new_offset); ehead_new = CacheCalloc(I->G,int,(map->Dim[0]*map->Dim[1]*map->Dim[2]), group_id,block_base + cCache_map_ehead_new_offset); #endif /* now do a filter-reassignment pass to remap fake vertices to the original line vertex while deleting duplicate entries */ memset( tempRef, 0, sizeof(int) * i_nVertex ); if(n_voxel < (3*n)) { /* faster to traverse the entire map */ int *start; for(a = iMin0; a <= iMax0; a++) { for(b = iMin1; b <= iMax1; b++) { for(c = iMin2; c <= iMax2; c++) { start = MapEStart(map,a,b,c); h = *start; if((h<0)&&(h>=neelem)) { sp = elist-h; *(start) = -h; /* flip sign */ i = *(sp++); if(ehead_new!=ehead) { ehead_new[(start - ehead)-1] = newelem; ip = ip0 = (elist_new+newelem); } else { ip = ip0 = (sp-1); } while(i>=0) { register int ii = *(sp++); if(i >= i_nVertex) /* reference -- remap */ i = tempRef[i]; if(!tempRef[i]) { /*eliminate duplicates */ tempRef[i] = 1; *(ip++) = i; } i = ii; } *(ip++) = -1; /* terminate list */ newelem+= (ip-ip0); /* now reset flags efficiently */ i = *(ip0++); while(i >= 0) { /* unroll */ ii = *(ip0++); tempRef[i] = 0; if((i=ii)>=0) { ii = *(ip0++); tempRef[i] = 0; if((i=ii)>=0) { ii = *(ip0++); tempRef[i] = 0; i=ii; } } } } } /* for c */ } /* for b */ } /* for a */ } else { int *start; int jm1,jp1,km1,kp1,lm1,lp1; MapType *map = I->Map; v = tempVertex; for(e = 0; e < n; e++) { MapLocus(map,v, &j, &k, &l); jm1 = j-1; jp1 = j+1; km1 = k-1; kp1 = k+1; lm1 = l-1; lp1 = l+1; if(perspective) { /* for normal maps */ register int *iPtr1 = map->EHead + ((j-1) * map->D1D2) + ((k-1) * map->Dim[2]) + (l-1); for(a = jm1; a <= jp1; a++) { register int *iPtr2 = iPtr1; if( (a >= iMin0) && (a <= iMax0) ) { for(b = km1; b <= kp1; b++) { register int *iPtr3 = iPtr2; if( (b >= iMin1) && (b <= iMax1) ) { for(c = lm1; c <= lp1; c++) { if( (c >= iMin2) && (c <= iMax2) ) { start = iPtr3; /*start = MapEStart(map,a,b,c);*/ h = *start; if((h<0)&&(h>=neelem)) { sp = elist - h; (*start) = -h; /* no repeat visits */ i = *(sp++); if(ehead_new!=ehead) { ehead_new[(start - ehead)] = newelem; ip = ip0 = (elist_new+newelem); } else { ip = ip0 = (sp-1); } while(i >= 0) { register int ii = *(sp++); if(i >= i_nVertex) i = tempRef[i]; if(!tempRef[i]) { /*eliminate duplicates */ tempRef[i] = 1; *(ip++) = i; } i = ii; } *(ip++) = -1; /* terminate list */ newelem+= (ip-ip0); /* now reset flags efficiently */ i = *(ip0++); while(i >= 0) { /* unroll */ ii = *(ip0++); tempRef[i] = 0; if((i=ii)>=0) { ii = *(ip0++); tempRef[i] = 0; if((i=ii)>=0) { ii = *(ip0++); tempRef[i] = 0; i=ii; } } } } /* h > 0 */ } iPtr3 ++; } } iPtr2 += map->Dim[2]; } } iPtr1 += map->D1D2; } } else { /* for XY maps... */ c = l; { register int *iPtr1 = ehead + ((j-1) * map->D1D2) + ((k-1) * map->Dim[2]) + c; if((c >= iMin2) && (c <= iMax2)) { for(a = jm1; a <= jp1; a++) { register int *iPtr2 = iPtr1; if( (a >= iMin0) && (a <= iMax0) ) { for(b = km1; b <= kp1; b++) { if( (b >= iMin1) && (b <= iMax1) ) { start = iPtr2; /*start = MapEStart(map,a,b,c);*/ h = *start; if((h<0)&&(h>=neelem)) { sp = elist - h; (*start) = -h; /* no repeat visits */ i = *(sp++); if(ehead_new!=ehead) { ehead_new[(start - ehead)] = newelem; ip = ip0 = (elist_new+newelem); } else { ip = ip0 = sp-1; } while(i >= 0) { register int ii = *(sp++); if(i >= i_nVertex) i = tempRef[i]; if(!tempRef[i]) { /*eliminate duplicates */ tempRef[i] = 1; *(ip++) = i; } i = ii; } *(ip++) = -1; /* terminate list */ newelem+= (ip-ip0); /* now reset flags efficiently */ i = *(ip0++); while(i >= 0) { /* unroll */ ii = *(ip0++); tempRef[i] = 0; if((i=ii)>=0) { ii = *(ip0++); tempRef[i] = 0; if((i=ii)>=0) { ii = *(ip0++); tempRef[i] = 0; i = ii; } } } } /* h > 0 */ } iPtr2 += map->Dim[2]; } /* for b */ } iPtr1 += map->D1D2; } /* for a */ } } } v += 3; /* happens for EVERY e! */ } /* for e */ } if(ehead!=ehead_new) { CacheFreeP(I->G,map->EHead,group_id,block_base + cCache_map_ehead_offset,false); VLACacheFreeP(I->G,map->EList,group_id,block_base + cCache_map_elist_offset,false); map->EList = elist_new; map->EHead = ehead_new; map->NEElem = newelem; MemoryCacheReplaceBlock(I->G, group_id, block_base + cCache_map_ehead_new_offset, block_base + cCache_map_ehead_offset); MemoryCacheReplaceBlock(I->G, group_id, block_base + cCache_map_elist_new_offset, block_base + cCache_map_elist_offset); VLACacheSize(G,map->EList,int,map->NEElem,group_id,block_base + cCache_map_elist_offset); } } CacheFreeP(I->G,tempVertex,group_id,cCache_basis_tempVertex,false); CacheFreeP(I->G,tempRef,group_id,cCache_basis_tempRef,false); } else { /* simple sphere mode */ I->Map = MapNewCached(I->G,-sep,I->Vertex,I->NVertex,NULL,group_id,block_base); if(perspective) { MapSetupExpressPerp(I->Map,I->Vertex,front,I->NVertex,false,NULL); } else { MapSetupExpressXYVert(I->Map,I->Vertex,I->NVertex,false); } } } /*========================================================================*/ void BasisInit(PyMOLGlobals *G,CBasis *I,int group_id) { I->G = G; I->Vertex = VLACacheAlloc(I->G,float,1,group_id,cCache_basis_vertex); I->Radius = VLACacheAlloc(I->G,float,1,group_id,cCache_basis_radius); I->Radius2 = VLACacheAlloc(I->G,float,1,group_id,cCache_basis_radius2); I->Normal = VLACacheAlloc(I->G,float,1,group_id,cCache_basis_normal); I->Vert2Normal = VLACacheAlloc(I->G,int,1,group_id,cCache_basis_vert2normal); I->Precomp = VLACacheAlloc(I->G,float,1,group_id,cCache_basis_precomp); I->Map=NULL; I->NVertex=0; I->NNormal=0; } /*========================================================================*/ void BasisFinish(CBasis *I,int group_id) { if(I->Map) { MapFree(I->Map); I->Map=NULL; } VLACacheFreeP(I->G,I->Radius2,group_id,cCache_basis_radius2,false); VLACacheFreeP(I->G,I->Radius,group_id,cCache_basis_radius,false); VLACacheFreeP(I->G,I->Vertex,group_id,cCache_basis_vertex,false); VLACacheFreeP(I->G,I->Vert2Normal,group_id,cCache_basis_vert2normal,false); VLACacheFreeP(I->G,I->Normal,group_id,cCache_basis_normal,false); VLACacheFreeP(I->G,I->Precomp,group_id,cCache_basis_precomp,false); I->Vertex=NULL; } /*========================================================================*/ void BasisTrianglePrecompute(float *v0,float *v1,float *v2,float *pre) { float det; subtract3f(v1,v0,pre); subtract3f(v2,v0,pre+3); det = pre[0]*pre[4] - pre[1]*pre[3]; if(fabs(det) < EPSILON) *(pre+6) = 0.0F; else { *(pre+6) = 1.0F; *(pre+7) = 1.0F/det; } } void BasisTrianglePrecomputePerspective(float *v0,float *v1,float *v2,float *pre) { subtract3f(v1,v0,pre); subtract3f(v2,v0,pre+3); } void BasisCylinderSausagePrecompute(float *dir,float *pre) { float ln = (float)(1.0f / sqrt1f(dir[1]*dir[1]+dir[0]*dir[0])); pre[0] = dir[1] * ln; pre[1] = -dir[0] * ln; } #else typedef int this_file_is_no_longer_empty; #endif