/* 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: -* -* -* Z* ------------------------------------------------------------------- */ #include"os_predef.h" #include"os_std.h" #include"os_gl.h" #include"Base.h" #include"MemoryDebug.h" #include"OOMac.h" #include"RepSurface.h" #include"Map.h" #include"Scene.h" #include"Sphere.h" #include"Setting.h" #include"Color.h" #include"ObjectMolecule.h" #include"Triangle.h" #include"Vector.h" #include"Feedback.h" #include"main.h" #include"Util.h" #include"CGO.h" #include"P.h" #include"Selector.h" #ifdef NT #undef NT #endif typedef struct RepSurface { Rep R; int N; int NT; int proximity; float *V,*VN,*VC; int *Vis; int *T,*S; /* S=strips */ int NDot; float *Dot; float *DotNormal; int solidFlag; int oneColorFlag,oneColor; int allVisibleFlag; int *LastVisib; int *LastColor; int Type; float max_vdw; CGO *debug; } RepSurface; void RepSurfaceRender(RepSurface *I,CRay *ray,Pickable **pick); void RepSurfaceFree(RepSurface *I); int RepSurfaceSameVis(RepSurface *I,CoordSet *cs); void RepSurfaceColor(RepSurface *I,CoordSet *cs); void RepSurfaceFree(RepSurface *I) { FreeP(I->V); FreeP(I->VN); FreeP(I->VC); FreeP(I->Vis); FreeP(I->LastColor); FreeP(I->LastVisib); CGOFree(I->debug); VLAFreeP(I->T); VLAFreeP(I->S); /* VLAFreeP(I->N);*/ OOFreeP(I); } void RepSurfaceGetSolventDots(RepSurface *I,CoordSet *cs, float probe_radius,SphereRec *sp, float *extent,int *present,int circumscribe); static int ZOrderFn(float *array,int l,int r) { return (array[l]<=array[r]); } static int ZRevOrderFn(float *array,int l,int r) { return (array[l]>=array[r]); } static int check_and_add(int *cache, int spacing, int t0,int t1) { int *rec; int cnt; t0++; t1++; rec = cache + spacing*t0; cnt=spacing; while(cnt>0) { if(*rec==t1) return 1; if(!*rec) { *rec = t1; break; } rec++; cnt--; } rec = cache + spacing*t1; cnt=spacing; while(cnt>0) { if(*rec==t0) return 1; if(!*rec) { *rec = t0; break; } rec++; cnt--; } return 0; } void RepSurfaceRender(RepSurface *I,CRay *ray,Pickable **pick) { float *v=I->V; float *vn=I->VN; float *vc=I->VC; int *t=I->T; int *s=I->S; int c=I->N; int *vi=I->Vis; float *col; float alpha; int t_mode; alpha = SettingGet_f(I->R.cs->Setting,I->R.obj->Setting,cSetting_transparency); alpha=1.0F-alpha; if(fabs(alpha-1.0)<R_SMALL4) alpha=1.0F; if(ray) { ray->fTransparentf(ray,1.0F-alpha); if(I->Type==1) { /* dot surface */ float radius; radius = SettingGet_f(I->R.cs->Setting,I->R.obj->Setting,cSetting_dot_radius); if(radius==0.0F) { radius = ray->PixelRadius*SettingGet_f(I->R.cs->Setting, I->R.obj->Setting, cSetting_dot_width)/1.4142F; } if(I->oneColorFlag) { ray->fColor3fv(ray,ColorGet(I->oneColor)); } if(c) while(c--) { if(*vi) { if(!I->oneColorFlag) { ray->fColor3fv(ray,vc); } ray->fSphere3fv(ray,v,radius); } vi++; vc+=3; v+=3; } } else if(I->Type==0) { /* solid surface */ c=I->NT; if(I->oneColorFlag) { col=ColorGet(I->oneColor); while(c--) { if((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2)))))) ray->fTriangle3fv(ray,v+(*t)*3,v+(*(t+1))*3,v+(*(t+2))*3, vn+(*t)*3,vn+(*(t+1))*3,vn+(*(t+2))*3, col,col,col); t+=3; } } else { while(c--) { if((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2)))))) if((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))) ray->fTriangle3fv(ray,v+(*t)*3,v+(*(t+1))*3,v+(*(t+2))*3, vn+(*t)*3,vn+(*(t+1))*3,vn+(*(t+2))*3, vc+(*t)*3,vc+(*(t+1))*3,vc+(*(t+2))*3); t+=3; } } } else if(I->Type==2) { /* triangle mesh surface */ float radius; int t0,t1,t2; int spacing = 10; int *cache = Calloc(int,spacing*(I->N+1)); radius = SettingGet_f(I->R.cs->Setting,I->R.obj->Setting,cSetting_mesh_radius); if(radius==0.0F) { radius = ray->PixelRadius*SettingGet_f(I->R.cs->Setting, I->R.obj->Setting, cSetting_mesh_width)/2.0F; } c=I->NT; if(I->oneColorFlag) { col=ColorGet(I->oneColor); while(c--) { t0 = (*t); t1 = (*(t+1)); t2 = (*(t+2)); if((I->proximity&&((*(vi+t0))||(*(vi+t1))||(*(vi+t2))))|| ((*(vi+t0))&&(*(vi+t1))&&(*(vi+t2)))) { if(!check_and_add(cache,spacing,t0,t1)) ray->fSausage3fv(ray,v+t0*3,v+t1*3,radius,col,col); if(!check_and_add(cache,spacing,t1,t2)) ray->fSausage3fv(ray,v+t1*3,v+t2*3,radius,col,col); if(!check_and_add(cache,spacing,t2,t0)) ray->fSausage3fv(ray,v+t2*3,v+t0*3,radius,col,col); } t+=3; } } else { while(c--) { t0 = (*t); t1 = (*(t+1)); t2 = (*(t+2)); if((I->proximity&&((*(vi+t0))||(*(vi+t1))||(*(vi+t2))))|| ((*(vi+t0))&&(*(vi+t1))&&(*(vi+t2)))) if((*(vi+t0))||(*(vi+t1))||(*(vi+t2))) { if(!check_and_add(cache,spacing,t0,t1)) ray->fSausage3fv(ray,v+t0*3,v+t1*3,radius,vc+t0*3,vc+t1*3); if(!check_and_add(cache,spacing,t1,t2)) ray->fSausage3fv(ray,v+t1*3,v+t2*3,radius,vc+t1*3,vc+t2*3); if(!check_and_add(cache,spacing,t2,t0)) ray->fSausage3fv(ray,v+t2*3,v+t0*3,radius,vc+t2*3,vc+t0*3); } t+=3; } } FreeP(cache); } ray->fTransparentf(ray,0.0); } else if(pick&&PMGUI) { } else if(PMGUI) { if(I->debug) CGORenderGL(I->debug,NULL,NULL,NULL); if(I->Type==1) { /* no triangle information, so we're rendering dots only */ int normals = SettingGet_f(I->R.cs->Setting,I->R.obj->Setting,cSetting_dot_normals); int lighting = SettingGet_f(I->R.cs->Setting,I->R.obj->Setting,cSetting_dot_lighting); int use_dlst; if(!normals) SceneResetNormal(true); if(!lighting) glDisable(GL_LIGHTING); use_dlst = (int)SettingGet(cSetting_use_display_lists); if(use_dlst&&I->R.displayList) { glCallList(I->R.displayList); } else { if(use_dlst) { if(!I->R.displayList) { I->R.displayList = glGenLists(1); if(I->R.displayList) { glNewList(I->R.displayList,GL_COMPILE_AND_EXECUTE); } } } glPointSize(SettingGet_f(I->R.cs->Setting,I->R.obj->Setting,cSetting_dot_width)); if(c) { glColor3f(1.0,0.0,0.0); glBegin(GL_POINTS); SceneResetNormal(true); if(I->oneColorFlag) { glColor3fv(ColorGet(I->oneColor)); } while(c--) { if(*vi) { if(!I->oneColorFlag) { glColor3fv(vc); } if(normals) glNormal3fv(vn); glVertex3fv(v); } vi++; vc+=3; vn+=3; v+=3; } glEnd(); } if(use_dlst&&I->R.displayList) { glEndList(); } if(!lighting) glEnable(GL_LIGHTING); } } else if(I->Type==2) { /* rendering triangle mesh */ int normals = SettingGet_b(I->R.cs->Setting,I->R.obj->Setting,cSetting_mesh_normals); int lighting = SettingGet_f(I->R.cs->Setting,I->R.obj->Setting,cSetting_mesh_lighting); int use_dlst; if(!normals) SceneResetNormal(true); if(!lighting) glDisable(GL_LIGHTING); use_dlst = (int)SettingGet(cSetting_use_display_lists); if(use_dlst&&I->R.displayList) { glCallList(I->R.displayList); } else { glLineWidth(SettingGet_f(I->R.cs->Setting,I->R.obj->Setting,cSetting_mesh_width)); if(use_dlst) { if(!I->R.displayList) { I->R.displayList = glGenLists(1); if(I->R.displayList) { glNewList(I->R.displayList,GL_COMPILE_AND_EXECUTE); } } } c=I->NT; if(c) { if(I->oneColorFlag) { glColor3fv(ColorGet(I->oneColor)); while(c--) { if((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2)))))) { if(normals) { glBegin(GL_LINE_STRIP); glNormal3fv(vn+(*(t+2))*3); glVertex3fv(v+(*(t+2))*3); glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glEnd(); } else { glBegin(GL_LINE_STRIP); glVertex3fv(v+(*(t+2))*3); glVertex3fv(v+(*t)*3); t++; glVertex3fv(v+(*t)*3); t++; glVertex3fv(v+(*t)*3); t++; glEnd(); } } else t+=3; } } else { while(c--) { if((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2)))))) { if(normals) { glBegin(GL_LINE_STRIP); glColor3fv(vc+(*(t+2))*3); glNormal3fv(vn+(*(t+2))*3); glVertex3fv(v+(*(t+2))*3); glColor3fv(vc+(*t)*3); glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glColor3fv(vc+(*t)*3); glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glColor3fv(vc+(*t)*3); glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glEnd(); } else { glBegin(GL_LINE_STRIP); glColor3fv(vc+(*(t+2))*3); glVertex3fv(v+(*(t+2))*3); glColor3fv(vc+(*t)*3); glVertex3fv(v+(*t)*3); t++; glColor3fv(vc+(*t)*3); glVertex3fv(v+(*t)*3); t++; glColor3fv(vc+(*t)*3); glVertex3fv(v+(*t)*3); t++; glEnd(); } } else t+=3; } } } if(use_dlst&&I->R.displayList) { glEndList(); } if(!lighting) glEnable(GL_LIGHTING); } } else { /* we're rendering triangles */ if(alpha!=1.0) { t_mode = SettingGet_i(I->R.cs->Setting,I->R.obj->Setting,cSetting_transparency_mode); if(t_mode) { float **t_buf=NULL,**tb; float *z_value=NULL,*zv; int *ix=NULL; int n_tri = 0; float sum[3]; float matrix[16]; glGetFloatv(GL_MODELVIEW_MATRIX,matrix); t_buf = Alloc(float*,I->NT*9); z_value = Alloc(float,I->NT); ix = Alloc(int,I->NT); zv = z_value; tb = t_buf; c = I->NT; if(I->oneColorFlag) { while(c--) { if((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2)))))) { *(tb++) = vn+(*t)*3; *(tb++) = v+(*t)*3; *(tb++) = vn+(*(t+1))*3; *(tb++) = v+(*(t+1))*3; *(tb++) = vn+(*(t+2))*3; *(tb++) = v+(*(t+2))*3; add3f(tb[-1],tb[-3],sum); add3f(sum,tb[-5],sum); *(zv++) = matrix[2]*sum[0]+matrix[6]*sum[1]+matrix[10]*sum[2]; n_tri++; } t+=3; } } else { while(c--) { if((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2)))))) if((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))) { *(tb++) = vc+(*t)*3; *(tb++) = vn+(*t)*3; *(tb++) = v+(*t)*3; *(tb++) = vc+(*(t+1))*3; *(tb++) = vn+(*(t+1))*3; *(tb++) = v+(*(t+1))*3; *(tb++) = vc+(*(t+2))*3; *(tb++) = vn+(*(t+2))*3; *(tb++) = v+(*(t+2))*3; add3f(tb[-1],tb[-4],sum); add3f(sum,tb[-7],sum); *(zv++) = matrix[2]*sum[0]+matrix[6]*sum[1]+matrix[10]*sum[2]; n_tri++; } t+=3; } } switch(t_mode) { case 1: UtilSortIndex(n_tri,z_value,ix,(UtilOrderFn*)ZOrderFn); break; default: UtilSortIndex(n_tri,z_value,ix,(UtilOrderFn*)ZRevOrderFn); break; } c=n_tri; if(I->oneColorFlag) { col=ColorGet(I->oneColor); glColor4f(col[0],col[1],col[2],alpha); glBegin(GL_TRIANGLES); for(c=0;c<n_tri;c++) { tb = t_buf+6*ix[c]; glNormal3fv(*(tb++)); glVertex3fv(*(tb++)); glNormal3fv(*(tb++)); glVertex3fv(*(tb++)); glNormal3fv(*(tb++)); glVertex3fv(*(tb++)); } glEnd(); } else { glBegin(GL_TRIANGLES); for(c=0;c<n_tri;c++) { float *vv; tb = t_buf+9*ix[c]; vv = *(tb++); glColor4f(vv[0],vv[1],vv[2],alpha); glNormal3fv(*(tb++)); glVertex3fv(*(tb++)); vv = *(tb++); glColor4f(vv[0],vv[1],vv[2],alpha); glNormal3fv(*(tb++)); glVertex3fv(*(tb++)); vv = *(tb++); glColor4f(vv[0],vv[1],vv[2],alpha); glNormal3fv(*(tb++)); glVertex3fv(*(tb++)); } glEnd(); } FreeP(ix); FreeP(z_value); FreeP(t_buf); } else { /* fast and ugly */ /* glCullFace(GL_BACK); glEnable(GL_CULL_FACE); glDepthMask(GL_FALSE);*/ if(I->allVisibleFlag) { if(I->oneColorFlag) { col = ColorGet(I->oneColor); glColor4f(col[0],col[1],col[2],alpha); c=*(s++); while(c) { glBegin(GL_TRIANGLE_STRIP); glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; while(c--) { glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; } glEnd(); c=*(s++); } } else { c=*(s++); while(c) { glBegin(GL_TRIANGLE_STRIP); col = vc+(*s)*3; glColor4f(col[0],col[1],col[2],alpha); glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; col = vc+(*s)*3; glColor4f(col[0],col[1],col[2],alpha); glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; while(c--) { col = vc+(*s)*3; glColor4f(col[0],col[1],col[2],alpha); glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; } glEnd(); c=*(s++); } } } else { /* subset s*/ c=I->NT; if(c) { glBegin(GL_TRIANGLES); if(I->oneColorFlag) { col = ColorGet(I->oneColor); glColor4f(col[0],col[1],col[2],alpha); while(c--) { if((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2)))))) { col = vc+(*t)*3; glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; col = vc+(*t)*3; glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; col = vc+(*t)*3; glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; } else t+=3; } } else { while(c--) { if((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2)))))) { col = vc+(*t)*3; glColor4f(col[0],col[1],col[2],alpha); glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; col = vc+(*t)*3; glColor4f(col[0],col[1],col[2],alpha); glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; col = vc+(*t)*3; glColor4f(col[0],col[1],col[2],alpha); glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; } else t+=3; } } glEnd(); } } /* glDisable(GL_CULL_FACE); glDepthMask(GL_TRUE);*/ } } else { /* opaque */ int use_dlst,simplify=0; use_dlst = (int)SettingGet(cSetting_use_display_lists); simplify = (int)SettingGet(cSetting_simplify_display_lists); if(use_dlst&&I->R.displayList) { glCallList(I->R.displayList); } else { if(use_dlst) { if(!I->R.displayList) { I->R.displayList = glGenLists(1); if(I->R.displayList) { glNewList(I->R.displayList,GL_COMPILE_AND_EXECUTE); } } } if(I->allVisibleFlag) { if(I->oneColorFlag) { if(use_dlst&&simplify) { /* simplify: try to help display list optimizer */ glColor3fv(ColorGet(I->oneColor)); c=*(s++); while(c) { glBegin(GL_TRIANGLES); s+=2; while(c--) { s-=2; glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; } glEnd(); c=*(s++); } } else { glColor3fv(ColorGet(I->oneColor)); c=*(s++); while(c) { glBegin(GL_TRIANGLE_STRIP); glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; while(c--) { glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; } glEnd(); c=*(s++); } } /* use_dlst&&simplify */ } else { /* not one color */ if(use_dlst&&simplify) { /* simplify: try to help display list optimizer */ c=*(s++); while(c) { glBegin(GL_TRIANGLES); s+=2; while(c--) { s-=2; glColor3fv(vc+(*s)*3); glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; glColor3fv(vc+(*s)*3); glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; glColor3fv(vc+(*s)*3); glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; } glEnd(); c=*(s++); } } else { c=*(s++); while(c) { glBegin(GL_TRIANGLE_STRIP); glColor3fv(vc+(*s)*3); glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; glColor3fv(vc+(*s)*3); glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; while(c--) { glColor3fv(vc+(*s)*3); glNormal3fv(vn+(*s)*3); glVertex3fv(v+(*s)*3); s++; } glEnd(); c=*(s++); } } } /* one color */ } else { /* subsets */ c=I->NT; if(c) { glBegin(GL_TRIANGLES); if(I->oneColorFlag) { glColor3fv(ColorGet(I->oneColor)); while(c--) { if((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2)))))) { glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; } else t+=3; } } else { while(c--) { if((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2)))))) { glColor3fv(vc+(*t)*3); glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glColor3fv(vc+(*t)*3); glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glColor3fv(vc+(*t)*3); glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; } else t+=3; } } glEnd(); } } if(use_dlst&&I->R.displayList) { glEndList(); } } } } if(SettingGet(cSetting_surface_debug)) { t=I->T; c=I->NT; if(c) { glBegin(GL_TRIANGLES); while(c--) { if(I->allVisibleFlag||((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2))))))) { glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; } else { t+=3; } } glEnd(); } t=I->T; c=I->NT; if(c) { glColor3f(0.0,1.0,0.0); while(c--) { glBegin(GL_LINE_STRIP); if(I->allVisibleFlag||((I->proximity&&((*(vi+(*t)))||(*(vi+(*(t+1))))||(*(vi+(*(t+2))))))|| ((*(vi+(*t)))&&(*(vi+(*(t+1))))&&(*(vi+(*(t+2))))))) { glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; glNormal3fv(vn+(*t)*3); glVertex3fv(v+(*t)*3); t++; } else { t+=3; } glEnd(); } } c=I->N; if(c) { glColor3f(1.0,0.0,0.0); glBegin(GL_LINES); SceneResetNormal(true); while(c--) { glVertex3fv(v); glVertex3f(v[0]+vn[0]/2,v[1]+vn[1]/2,v[2]+vn[2]/2); v+=3; vn+=3; } glEnd(); } glColor3f(0.3F,0.3F,1.0F); v=TestLine; c=NTestLine; glBegin(GL_LINES); while(c--) { glVertex3fv(v); v+=3; glVertex3fv(v); v+=3; } glEnd(); } glColor4f(1,1,1,1); SceneResetNormal(0); v=TestLine; c=NTestLine; glBegin(GL_LINES); while(c--) { glVertex3fv(v); v+=3; glVertex3fv(v); v+=3; } glEnd(); } } int RepSurfaceSameVis(RepSurface *I,CoordSet *cs) { int same = true; int *lv,*lc,*cc; int a; AtomInfoType *ai; ai = cs->Obj->AtomInfo; lv = I->LastVisib; lc = I->LastColor; cc = cs->Color; for(a=0;a<cs->NIndex;a++) { if(*(lv++)!=(ai + cs->IdxToAtm[a])->visRep[cRepSurface] ) { same=false; break; } if(*(lc++)!=*(cc++)) { same=false; break; } } return(same); } void RepSurfaceColor(RepSurface *I,CoordSet *cs) { MapType *map; int a,i0,i,j,c1; float *v0,*vc,*c0; float *n0; int *vi,*lv,*lc,*cc; int first_color; float *v_pos,v_above[3]; int ramp_above; ObjectMolecule *obj; float probe_radius; float dist,minDist; float cutoff; int inclH; int cullByFlag = false; int surface_mode; int surface_color; int *present=NULL,*ap; int carve_state = 0; int carve_flag = false; float carve_cutoff; char *carve_selection = NULL; float *carve_vla = NULL; MapType *carve_map = NULL; int clear_state = 0; int clear_flag = false; float clear_cutoff; char *clear_selection = NULL; float *clear_vla = NULL; MapType *clear_map = NULL; AtomInfoType *ai2,*ai1; obj=cs->Obj; surface_mode = SettingGet_i(cs->Setting,obj->Obj.Setting,cSetting_surface_mode); ramp_above = SettingGet_i(cs->Setting,obj->Obj.Setting,cSetting_surface_ramp_above_mode); surface_color = SettingGet_color(cs->Setting,obj->Obj.Setting,cSetting_surface_color); cullByFlag = (surface_mode==cRepSurface_by_flags); inclH = !(surface_mode==cRepSurface_heavy_atoms); probe_radius = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_solvent_radius); I->proximity = SettingGet_b(cs->Setting,obj->Obj.Setting,cSetting_surface_proximity); carve_cutoff = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_carve_cutoff); clear_cutoff = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_clear_cutoff); cutoff = I->max_vdw+2*probe_radius; if(!I->LastVisib) I->LastVisib = Alloc(int,cs->NIndex); if(!I->LastColor) I->LastColor = Alloc(int,cs->NIndex); lv = I->LastVisib; lc = I->LastColor; cc = cs->Color; ai2=obj->AtomInfo; for(a=0;a<cs->NIndex;a++) { *(lv++) = (ai2 + cs->IdxToAtm[a])->visRep[cRepSurface]; *(lc++) = *(cc++); } if(I->N) { if(carve_cutoff>0.0F) { carve_state = SettingGet_i(cs->Setting,obj->Obj.Setting,cSetting_surface_carve_state) - 1; carve_selection = SettingGet_s(cs->Setting,obj->Obj.Setting,cSetting_surface_carve_selection); if(carve_selection) carve_map = SelectorGetSpacialMapFromSeleCoord(SelectorIndexByName(carve_selection), carve_state, carve_cutoff,&carve_vla); if(carve_map) MapSetupExpress(carve_map); carve_flag = true; } if(clear_cutoff>0.0F) { clear_state = SettingGet_i(cs->Setting,obj->Obj.Setting,cSetting_surface_clear_state) - 1; clear_selection = SettingGet_s(cs->Setting,obj->Obj.Setting,cSetting_surface_clear_selection); if(clear_selection) clear_map = SelectorGetSpacialMapFromSeleCoord(SelectorIndexByName(clear_selection), clear_state, clear_cutoff,&clear_vla); if(clear_map) MapSetupExpress(clear_map); clear_flag = true; } if(!I->VC) I->VC = Alloc(float,3*I->N); vc=I->VC; if(!I->Vis) I->Vis = Alloc(int,I->N); vi=I->Vis; if(ColorCheckRamped(surface_color)) { I->oneColorFlag=false; } else { I->oneColorFlag=true; } first_color=-1; present = Alloc(int,cs->NIndex); ap = present; for(a=0;a<cs->NIndex;a++) { ai1 = obj->AtomInfo + cs->IdxToAtm[a]; if(ai1->visRep[cRepSurface]&& (inclH||(!ai1->hydrogen))&& ((!cullByFlag)|| (!(ai1->flags&(cAtomFlag_ignore|cAtomFlag_exfoliate))))) *ap = 2; else *ap = 0; ap++; } map=MapNewFlagged(2*I->max_vdw+probe_radius,cs->Coord,cs->NIndex,NULL,present); MapSetupExpress(map); for(a=0;a<cs->NIndex;a++) if(!present[a]) { ai1 = obj->AtomInfo+cs->IdxToAtm[a]; if((!cullByFlag)||!(ai1->flags&cAtomFlag_ignore)) { v0 = cs->Coord+3*a; i=*(MapLocusEStart(map,v0)); if(i) { j=map->EList[i++]; while(j>=0) { if(present[j]>1) { ai2 = obj->AtomInfo+cs->IdxToAtm[j]; if(within3f(cs->Coord+3*j,v0,ai1->vdw+ai2->vdw+probe_radius)) { present[a]=1; break; } } j=map->EList[i++]; } } } } MapFree(map); map = NULL; /* now, assign colors to each point */ map=MapNewFlagged(cutoff,cs->Coord,cs->NIndex,NULL,present); if(map) { MapSetupExpress(map); for(a=0;a<I->N;a++) { c1=1; minDist=MAXFLOAT; i0=-1; v0 = I->V+3*a; n0 = I->VN+3*a; vi = I->Vis+a; /* colors */ i=*(MapLocusEStart(map,v0)); if(i) { j=map->EList[i++]; while(j>=0) { ai2 = obj->AtomInfo + cs->IdxToAtm[j]; if((inclH||(!ai2->hydrogen))&& ((!cullByFlag)|| (!(ai2->flags&cAtomFlag_ignore)))) { dist = (float)diff3f(v0,cs->Coord+j*3)-ai2->vdw; if(dist<minDist) { i0=j; minDist=dist; } } j=map->EList[i++]; } } if(i0>=0) { c1=*(cs->Color+i0); if(I->oneColorFlag) { if(first_color>=0) { if(first_color!=c1) I->oneColorFlag=false; } else first_color=c1; } if(I->allVisibleFlag) *vi = 1; else { ai2 = obj->AtomInfo+cs->IdxToAtm[i0]; if(ai2->visRep[cRepSurface]&& (inclH||(!ai2->hydrogen))&& ((!cullByFlag)|| (!(ai2->flags&(cAtomFlag_ignore|cAtomFlag_exfoliate))))) *vi = 1; else *vi = 0; } } else { *vi = 0; } if(carve_flag && (*vi)) { /* is point visible, and are we carving? */ *vi = 0; if(carve_map) { minDist=MAXFLOAT; i=*(MapLocusEStart(carve_map,v0)); if(i) { j=carve_map->EList[i++]; while(j>=0) { if(within3f(carve_vla+3*j,v0,carve_cutoff)) { *vi=1; break; } j=carve_map->EList[i++]; } } } } if(clear_flag && (*vi)) { /* is point visible, and are we clearing? */ if(clear_map) { minDist=MAXFLOAT; i=*(MapLocusEStart(clear_map,v0)); if(i) { j=clear_map->EList[i++]; while(j>=0) { if(within3f(clear_vla+3*j,v0,clear_cutoff)) { *vi=0; break; } j=clear_map->EList[i++]; } } } } if(ColorCheckRamped(surface_color)) { c1 = surface_color; } if(ColorCheckRamped(c1)) { I->oneColorFlag=false; switch(ramp_above) { case 1: copy3f(n0,v_above); scale3f(v_above,probe_radius,v_above); add3f(v0,v_above,v_above); v_pos = v_above; break; default: v_pos = v0; break; } ColorGetRamped(c1,v_pos,vc); vc+=3; } else { c0 = ColorGet(c1); *(vc++) = *(c0++); *(vc++) = *(c0++); *(vc++) = *(c0++); } if(!*vi) I->allVisibleFlag=false; vi++; } MapFree(map); } if(I->oneColorFlag) { I->oneColor=first_color; } } if(surface_color>=0) { I->oneColorFlag=true; I->oneColor=surface_color; } if(PMGUI) { if(I->R.displayList) { if(PIsGlutThread()) { glDeleteLists(I->R.displayList,1); I->R.displayList = 0; } else { char buffer[255]; /* pass this off to the main thread */ sprintf(buffer,"_cmd.gl_delete_lists(%d,%d)\n",I->R.displayList,1); PParse(buffer); } } } if(carve_map) MapFree(carve_map); VLAFreeP(carve_vla); if(clear_map) MapFree(clear_map); VLAFreeP(clear_vla); FreeP(present); } Rep *RepSurfaceNew(CoordSet *cs) { ObjectMolecule *obj; int a,b,i,j,c; MapType *map,*solv_map; float *v0=NULL,*v,*vn=NULL,*vn0=NULL,*extent=NULL,*n0=NULL; int SurfaceFlag = false; float probe_radius,probe_radius2; float probe_rad_more,probe_rad_more2; float probe_rad_less,probe_rad_less2; float trim_cutoff,trim_factor; int inclH = true; int cullByFlag = false; int flag,*dot_flag,*p; float minimum_sep; int visFlag; int surface_quality; int surface_mode; int *present = NULL,*ap; int pres_flag; int surface_type; int surface_solvent; int MaxN; SphereRec *sp = Sphere0; SphereRec *ssp = Sphere0; AtomInfoType *ai1,*ai2; int n_present = 0; float solv_tole; int carve_state = 0; int carve_flag = false; float carve_cutoff; char *carve_selection = NULL; float *carve_vla = NULL; MapType *carve_map = NULL; int circumscribe = 0; #if 0 int c1; float v1[3]; float vdw; #endif OOAlloc(RepSurface); obj = cs->Obj; I->max_vdw = ObjectMoleculeGetMaxVDW(obj); surface_mode = SettingGet_i(cs->Setting,obj->Obj.Setting,cSetting_surface_mode); surface_type = SettingGet_i(cs->Setting,obj->Obj.Setting,cSetting_surface_type); surface_solvent = SettingGet_b(cs->Setting,obj->Obj.Setting,cSetting_surface_solvent); I->Type = surface_type; trim_cutoff = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_trim_cutoff); trim_factor = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_trim_factor); cullByFlag = (surface_mode==cRepSurface_by_flags); inclH = !(surface_mode==cRepSurface_heavy_atoms); visFlag=false; for(a=0;a<cs->NIndex;a++) { ai1=obj->AtomInfo+cs->IdxToAtm[a]; if(ai1->visRep[cRepSurface]&& (inclH||(!ai1->hydrogen))&& ((!cullByFlag)| (!(ai1->flags&(cAtomFlag_exfoliate|cAtomFlag_ignore))))) { visFlag=true; break; } } if(!visFlag) { OOFreeP(I); return(NULL); /* skip if no thing visible */ } RepInit(&I->R); surface_quality = SettingGet_i(cs->Setting,obj->Obj.Setting,cSetting_surface_quality); if(surface_quality>=4) { /* totally impractical */ minimum_sep = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_best)/4; sp=Sphere4; ssp=Sphere4; circumscribe = 120; } else if(surface_quality>=3) { /* nearly impractical */ minimum_sep = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_best)/3; sp=Sphere4; ssp=Sphere3; circumscribe = 90; } else if(surface_quality>=2) { /* nearly perfect */ minimum_sep = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_best)/2; sp=Sphere3; ssp=Sphere3; circumscribe = 60; } else if(surface_quality>=1) { /* good */ minimum_sep = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_best); sp=Sphere2; ssp=Sphere3; } else if(!surface_quality) { /* 0 - normal */ minimum_sep = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_normal); sp=Sphere1; ssp=Sphere2; } else if(surface_quality==-1) { /* -1 */ minimum_sep = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_poor); sp=Sphere1; ssp=Sphere2; } else if(surface_quality==-2) { /* -2 god awful*/ minimum_sep = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_poor)*1.5F; sp=Sphere1; ssp=Sphere1; } else if(surface_quality==-3) { /* -3 miserable */ minimum_sep = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_miserable); sp=Sphere1; ssp=Sphere1; } else { minimum_sep = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_miserable)*1.18F; sp=Sphere0; ssp=Sphere1; } probe_radius = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_solvent_radius); if(!surface_solvent) { if(probe_radius<(2.5F*minimum_sep)) { probe_radius = 2.5F*minimum_sep; } } else { circumscribe = 0; } solv_tole = minimum_sep * 0.04; probe_radius2 = probe_radius*probe_radius; probe_rad_more = probe_radius*(1.0F+solv_tole); probe_rad_more2 = probe_rad_more * probe_rad_more; if(surface_type!=0) { /* not a solid surface */ probe_rad_less = probe_radius*(1.0F-solv_tole); } else { /* solid surface */ if(surface_quality>2) { probe_rad_less = probe_radius*(1.0F-solv_tole/2); } else if(surface_quality>1) { probe_rad_less = probe_radius*(1.0F-solv_tole); } else { probe_rad_less = probe_radius; } } probe_rad_less2 = probe_rad_less * probe_rad_less; I->N=0; I->NT=0; I->S=NULL; I->V=NULL; I->VC=NULL; I->Vis=NULL; I->VN=NULL; I->T=NULL; I->Dot=NULL; I->NDot=0; I->LastVisib=NULL; I->LastColor=NULL; I->R.fRender=(void (*)(struct Rep *, CRay *, Pickable **))RepSurfaceRender; I->R.fFree=(void (*)(struct Rep *))RepSurfaceFree; I->R.fRecolor=(void (*)(struct Rep*, struct CoordSet*))RepSurfaceColor; I->R.fSameVis=(int (*)(struct Rep*, struct CoordSet*))RepSurfaceSameVis; I->R.obj = (CObject*)(cs->Obj); I->R.cs = cs; I->allVisibleFlag=true; I->debug = NULL; obj = cs->Obj; /* don't waist time computing a Surface unless we need it!! */ for(a=0;a<cs->NIndex;a++) { ai1=obj->AtomInfo+cs->IdxToAtm[a]; if(ai1->visRep[cRepSurface]&& ((!cullByFlag)| (!(ai1->flags&(cAtomFlag_exfoliate))))) SurfaceFlag=true; else I->allVisibleFlag=false; } if(SurfaceFlag) { OrthoBusyFast(0,1); n_present = cs->NIndex; carve_selection = SettingGet_s(cs->Setting,obj->Obj.Setting,cSetting_surface_carve_selection); carve_cutoff = SettingGet_f(cs->Setting,obj->Obj.Setting,cSetting_surface_carve_cutoff); if((!carve_selection)||(!carve_selection[0])) carve_cutoff=0.0F; if(carve_cutoff>0.0F) { carve_state = SettingGet_i(cs->Setting,obj->Obj.Setting,cSetting_surface_carve_state) - 1; carve_cutoff += 2*I->max_vdw+probe_radius; if(carve_selection) carve_map = SelectorGetSpacialMapFromSeleCoord(SelectorIndexByName(carve_selection), carve_state, carve_cutoff, &carve_vla); if(carve_map) MapSetupExpress(carve_map); carve_flag = true; I->allVisibleFlag=false; } if(!I->allVisibleFlag) { /* optimize the space over which we calculate a surface */ /* first find out which atoms are really present */ present = Alloc(int,cs->NIndex); ap = present; for(a=0;a<cs->NIndex;a++) { ai1 = obj->AtomInfo + cs->IdxToAtm[a]; if(ai1->visRep[cRepSurface]&& (inclH||(!ai1->hydrogen))&& ((!cullByFlag)|| (!(ai1->flags&(cAtomFlag_ignore|cAtomFlag_exfoliate))))) *ap = 2; else *ap = 0; ap++; } map=MapNewFlagged(2*I->max_vdw+probe_radius,cs->Coord,cs->NIndex,extent,present); MapSetupExpress(map); for(a=0;a<cs->NIndex;a++) if(!present[a]) { ai1 = obj->AtomInfo+cs->IdxToAtm[a]; if((!cullByFlag)||!(ai1->flags&cAtomFlag_ignore)) { v0 = cs->Coord+3*a; i=*(MapLocusEStart(map,v0)); if(i) { j=map->EList[i++]; while(j>=0) { if(present[j]>1) { ai2 = obj->AtomInfo+cs->IdxToAtm[j]; if(within3f(cs->Coord+3*j,v0,ai1->vdw+ai2->vdw+probe_radius)) { present[a]=1; break; } } j=map->EList[i++]; } } } } if(carve_flag) { for(a=0;a<cs->NIndex;a++) { int include_flag = false; if(carve_map) { v0 = cs->Coord+3*a; i=*(MapLocusEStart(carve_map,v0)); if(i) { j=carve_map->EList[i++]; while(j>=0) { if(within3f(carve_vla+3*j,v0,carve_cutoff)) { include_flag = true; break; } j=carve_map->EList[i++]; } } } if(!include_flag) present[a]=0; } } MapFree(map); map = NULL; /* now count how many atoms we actually need to think about */ n_present = 0; for(a=0;a<cs->NIndex;a++) if(present[a]) { n_present++; } } if(n_present<1) n_present=1; /* safety */ MaxN = n_present*sp->nDot*10; I->V=Alloc(float,(MaxN+1)*3); ErrChkPtr(I->V); I->VN=Alloc(float,(MaxN+1)*3); ErrChkPtr(I->VN); I->N=0; v=I->V; vn=I->VN; RepSurfaceGetSolventDots(I,cs,probe_radius,ssp,extent,present,circumscribe); if(!surface_solvent) { map=MapNewFlagged(I->max_vdw+probe_rad_more,cs->Coord,cs->NIndex,extent,present); solv_map=MapNew(probe_rad_less,I->Dot,I->NDot,extent); /* I->debug=CGONew(); CGOBegin(I->debug,GL_POINTS); for(a=0;a<I->NDot;a++) CGOVertexv(I->debug,I->Dot+3*a); CGOEnd(I->debug);*/ if(map&&solv_map) { MapSetupExpress(solv_map); MapSetupExpress(map); if(I->NDot) { Vector3f *dot = NULL; dot=Alloc(Vector3f,sp->nDot); for(b=0;b<sp->nDot;b++) { scale3f(sp->dot[b],probe_radius,dot[b]); } v0 = I->Dot; for(a=0;a<I->NDot;a++) { OrthoBusyFast(a+I->NDot*2,I->NDot*5); /* 2/5 to 3/5 */ for(b=0;b<sp->nDot;b++) { register int ii; v[0]=v0[0]+dot[b][0]; v[1]=v0[1]+dot[b][1]; v[2]=v0[2]+dot[b][2]; flag=true; ii=*(MapLocusEStart(solv_map,v)); if(ii) { register int jj; register int *elist = solv_map->EList; register float *i_dot = I->Dot; register float v_0=v[0], v_1=v[1], v_2=v[2]; register float dist=probe_rad_less; register float dist2=probe_rad_less2; register float *v1,dx,dy,dz; jj=elist[ii++]; v1 = i_dot + 3*jj; while(jj>=0) { if(jj!=a) { /* huge bottleneck -- optimized for superscaler processors */ dx = v1[0]-v_0; dy = v1[1]-v_1; dz = v1[2]-v_2; dx = fabs(dx); dy = fabs(dy); if(!(dx>dist)) { dx = dx * dx; if(!(dy>dist)) { dz = fabs(dz); dy = dy * dy; if(!(dz>dist)) { dx = dx + dy; dz = dz * dz; if(!(dx>dist2)) if((dx + dz)<=dist2) { flag = false; break; } } } } } jj=elist[ii++]; v1 = i_dot + 3*jj; } } if(flag) { i=*(MapLocusEStart(map,v)); if(i) { j=map->EList[i++]; while(j>=0) { ai2 = obj->AtomInfo + cs->IdxToAtm[j]; if(present) pres_flag=present[j]; else pres_flag=((inclH||(!ai2->hydrogen))&& ((!cullByFlag)|| (!(ai2->flags&cAtomFlag_ignore)))); if(pres_flag) if(within3f(cs->Coord+3*j,v,ai2->vdw+probe_rad_more)) { flag=false; break; } j=map->EList[i++]; } } if(!flag) { vn[0]=-sp->dot[b][0]; vn[1]=-sp->dot[b][1]; vn[2]=-sp->dot[b][2]; if(I->N<MaxN) { I->N++; v+=3; vn+=3; } } } } v0 +=3; } FreeP(dot); } MapFree(solv_map); MapFree(map); } } else { v0 = I->Dot; n0 = I->DotNormal; if(I->NDot) { for(a=0;a<I->NDot;a++) { *(v++)=*(v0++); *(vn++)=*(n0++); *(v++)=*(v0++); *(vn++)=*(n0++); *(v++)=*(v0++); *(vn++)=*(n0++); I->N++; } } } FreeP(I->Dot); FreeP(I->DotNormal); /* now, eliminate dots that are too close to each other*/ /* CGOColor(I->debug,0.0,1.0,0.0); CGOBegin(I->debug,GL_POINTS); for(a=0;a<I->N;a++) CGOVertexv(I->debug,I->V+3*a); CGOEnd(I->debug); */ PRINTFB(FB_RepSurface,FB_Details) " RepSurface: %i surface points.\n",I->N ENDFB; if(I->N) { int repeat_flag=true; dot_flag=Alloc(int,I->N); while(repeat_flag) { repeat_flag=false; for(a=0;a<I->N;a++) dot_flag[a]=1; map=MapNew(minimum_sep,I->V,I->N,extent); MapSetupExpress(map); v=I->V; vn=I->VN; for(a=0;a<I->N;a++) { if(dot_flag[a]) { i=*(MapLocusEStart(map,v)); if(i) { j=map->EList[i++]; while(j>=0) { if(j!=a) { if(dot_flag[j]) { if(within3f(I->V+(3*j),v,minimum_sep)) { dot_flag[j]=0; add3f(vn,I->VN+(3*j),vn); average3f(I->V+(3*j),v,v); repeat_flag=true; } } } j=map->EList[i++]; } } } v+=3; vn+=3; } MapFree(map); v0=I->V; v=I->V; vn0=I->VN; vn=I->VN; p=dot_flag; c=I->N; I->N=0; for(a=0;a<c;a++) { if(*(p++)) { *(v0++)=*(v++); *(v0++)=*(v++); *(v0++)=*(v++); normalize3f(vn); *(vn0++)=*(vn++); *(vn0++)=*(vn++); *(vn0++)=*(vn++); I->N++; } else { v+=3; vn+=3; } } } FreeP(dot_flag); if(I->N) { I->V = ReallocForSure(I->V,float,(v0-I->V)); I->VN = ReallocForSure(I->VN,float,(vn0-I->VN)); } } /* now eliminate troublesome vertices in regions of extremely high curvature */ if(I->N&&(trim_cutoff>0.0F)&&(trim_factor>0.0F)) { int repeat_flag=true; float neighborhood = trim_factor*minimum_sep; float *v0,dot_sum; int n_nbr; dot_flag=Alloc(int,I->N); while(repeat_flag) { repeat_flag=false; for(a=0;a<I->N;a++) dot_flag[a]=1; map=MapNew(neighborhood,I->V,I->N,extent); MapSetupExpress(map); v=I->V; vn=I->VN; for(a=0;a<I->N;a++) { if(dot_flag[a]) { i=*(MapLocusEStart(map,v)); if(i) { n_nbr = 0; dot_sum = 0.0F; j=map->EList[i++]; while(j>=0) { if(j!=a) { if(dot_flag[j]) { v0 = I->V+3*j; if(within3f(v0,v,neighborhood)) { n0 = I->VN + 3*j; dot_sum += dot_product3f(n0,vn); n_nbr++; } } } j=map->EList[i++]; } if(n_nbr) { dot_sum/=n_nbr; if(dot_sum<trim_cutoff) { dot_flag[a]=false; repeat_flag = true; } } } } v+=3; vn+=3; } MapFree(map); v0=I->V; v=I->V; vn0=I->VN; vn=I->VN; p=dot_flag; c=I->N; I->N=0; for(a=0;a<c;a++) { if(*(p++)) { *(v0++)=*(v++); *(v0++)=*(v++); *(v0++)=*(v++); normalize3f(vn); *(vn0++)=*(vn++); *(vn0++)=*(vn++); *(vn0++)=*(vn++); I->N++; } else { v+=3; vn+=3; } } } FreeP(dot_flag); if(I->N) { I->V = ReallocForSure(I->V,float,(v0-I->V)); I->VN = ReallocForSure(I->VN,float,(vn0-I->VN)); } } PRINTFD(FB_RepSurface) " RepSurfaceNew-DEBUG: %i surface points after trimming.\n",I->N ENDFD; RepSurfaceColor(I,cs); PRINTFD(FB_RepSurface) " RepSurfaceNew-DEBUG: %i surface points after coloring.\n",I->N ENDFD; OrthoBusyFast(3,5); if(I->N) { if(surface_type!=1) { /* not a dot surface... */ float cutoff = minimum_sep*5.0F; if((cutoff>probe_radius)&&(!surface_solvent)) cutoff = probe_radius; I->T=TrianglePointsToSurface(I->V,I->VN,I->N,cutoff,&I->NT,&I->S,extent); PRINTFB(FB_RepSurface,FB_Details) " RepSurface: %i triangles.\n",I->NT ENDFB; } } else { I->V = ReallocForSure(I->V,float,1); I->VN = ReallocForSure(I->VN,float,1); } } if(carve_map) MapFree(carve_map); VLAFreeP(carve_vla); if(I->debug) CGOStop(I->debug); OrthoBusyFast(4,4); FreeP(present); return((void*)(struct Rep*)I); } void RepSurfaceGetSolventDots(RepSurface *I,CoordSet *cs, float probe_radius,SphereRec *sp, float *extent,int *present, int circumscribe) { ObjectMolecule *obj; int a,b,c=0,flag,i,j,ii,jj; float *v,*v0,vdw,*v1,*v2; float *n,*n0; MapType *map=NULL,*map2=NULL; int *p,*dot_flag; int cavity_cull,skip_flag; float probe_radius_plus; int dotCnt=0,maxCnt,maxDot=0; int cnt; int surface_mode; int surface_solvent; int cullByFlag; int inclH; int pres_flag; int stopDot; AtomInfoType *ai1,*ai2,*ai3; obj = cs->Obj; surface_mode = SettingGet_i(cs->Setting,obj->Obj.Setting,cSetting_surface_mode); cullByFlag = (surface_mode==cRepSurface_by_flags); inclH = !(surface_mode==cRepSurface_heavy_atoms); surface_solvent = SettingGet_b(cs->Setting,obj->Obj.Setting,cSetting_surface_solvent); cavity_cull = SettingGet_i(cs->Setting,obj->Obj.Setting,cSetting_cavity_cull); stopDot = cs->NIndex*sp->nDot+2*circumscribe; I->Dot=(float*)mmalloc(sizeof(float)*(stopDot+1)*3); I->DotNormal=(float*)mmalloc(sizeof(float)*(stopDot+1)*3); ErrChkPtr(I->Dot); ErrChkPtr(I->DotNormal); probe_radius_plus = probe_radius * 1.5F; I->NDot=0; map=MapNewFlagged(I->max_vdw+probe_radius,cs->Coord,cs->NIndex,extent,present); if(map) { MapSetupExpress(map); maxCnt=0; v=I->Dot; n=I->DotNormal; for(a=0;a<cs->NIndex;a++) { OrthoBusyFast(a,cs->NIndex*5); ai1 = obj->AtomInfo+cs->IdxToAtm[a]; if(present) pres_flag=present[a]; else pres_flag = (inclH||(!ai1->hydrogen))&& ((!cullByFlag)|| (!(ai1->flags&(cAtomFlag_ignore)))); if(pres_flag) { dotCnt=0; v0 = cs->Coord+3*a; vdw = ai1->vdw+probe_radius; skip_flag=false; i=*(MapLocusEStart(map,v0)); if(i) { j=map->EList[i++]; while(j>=0) { ai2 = obj->AtomInfo + cs->IdxToAtm[j]; if(j>a) /* only check if this is atom trails */ if((inclH||(!ai2->hydrogen))&& ((!cullByFlag)|| (!(ai2->flags&cAtomFlag_ignore)))) { if((ai2->vdw == ai1->vdw)) { /* handle singularities */ v1 = cs->Coord+3*j; if((v0[0]==v1[0]) && (v0[1]==v1[1]) && (v0[2]==v1[2])) skip_flag=true; } } j=map->EList[i++]; } } if(!skip_flag) { for(b=0;b<sp->nDot;b++) { v[0]=v0[0]+vdw*(n[0] = sp->dot[b][0]); v[1]=v0[1]+vdw*(n[1] = sp->dot[b][1]); v[2]=v0[2]+vdw*(n[2] = sp->dot[b][2]); flag=true; i=*(MapLocusEStart(map,v)); if(i) { j=map->EList[i++]; while(j>=0) { ai2 = obj->AtomInfo + cs->IdxToAtm[j]; if((inclH||(!ai2->hydrogen))&& ((!cullByFlag)|| (!(ai2->flags&cAtomFlag_ignore)))) if(j!=a) { skip_flag=false; if(ai1->vdw==ai2->vdw) { /* handle singularities */ v1 = cs->Coord+3*j; if((v0[0]==v1[0]) && (v0[1]==v1[1]) && (v0[2]==v1[2])) skip_flag=true; } if(!skip_flag) if(within3f(cs->Coord+3*j,v,ai2->vdw+probe_radius)) { flag=false; break; } } j=map->EList[i++]; } } if(flag&&(dotCnt<stopDot)) { dotCnt++; v+=3; n+=3; I->NDot++; } } } if(dotCnt>maxCnt) { maxCnt=dotCnt; maxDot=I->NDot-1; } } } /* for each pair of proximal atoms, circumscribe a circle for their intersection */ if(circumscribe&&(!surface_solvent)) map2=MapNewFlagged(2*(I->max_vdw+probe_radius),cs->Coord,cs->NIndex,extent,present); if(map2) { MapSetupExpress(map2); for(a=0;a<cs->NIndex;a++) { ai1 = obj->AtomInfo+cs->IdxToAtm[a]; if(present) pres_flag=present[a]; else pres_flag = (inclH||(!ai1->hydrogen))&& ((!cullByFlag)|| (!(ai1->flags&(cAtomFlag_ignore)))); if(pres_flag) { float vdw2; v0 = cs->Coord+3*a; vdw = ai1->vdw+probe_radius; vdw2 = vdw*vdw; skip_flag=false; i=*(MapLocusEStart(map2,v0)); if(i) { j=map2->EList[i++]; while(j>=0) { ai3 = obj->AtomInfo + cs->IdxToAtm[j]; if(j>a) /* only check if this is atom trails */ if((inclH||(!ai3->hydrogen))&& ((!cullByFlag)|| (!(ai3->flags&cAtomFlag_ignore)))) { if((ai3->vdw == ai1->vdw)) { /* handle singularities */ v2 = cs->Coord+3*j; if((v0[0]==v2[0]) && (v0[1]==v2[1]) && (v0[2]==v2[2])) skip_flag=true; } } j=map2->EList[i++]; } } if(!skip_flag) { ii=*(MapLocusEStart(map2,v0)); if(ii) { jj=map2->EList[ii++]; while(jj>=0) { float dist; ai3 = obj->AtomInfo + cs->IdxToAtm[jj]; if(jj>a) /* only check if this is atom trails */ if((inclH||(!ai3->hydrogen))&& ((!cullByFlag)|| (!(ai3->flags&cAtomFlag_ignore)))) { float vdw3 = ai3->vdw+probe_radius; v2 = cs->Coord+3*jj; dist = diff3f(v0,v2); if((dist>R_SMALL4)&&(dist<(vdw+vdw3))) { float vz[3],vx[3],vy[3], vp[3]; float tri_a=vdw, tri_b=vdw3, tri_c = dist; float tri_s = (tri_a+tri_b+tri_c)*0.5F; float area = sqrt1f(tri_s*(tri_s-tri_a)* (tri_s-tri_b)*(tri_s-tri_c)); float radius = (2*area)/dist; float adj = sqrt1f(vdw2 - radius*radius); subtract3f(v2,v0,vz); get_system1f3f(vz,vx,vy); copy3f(vz,vp); scale3f(vp,adj,vp); add3f(v0,vp,vp); for(b=0;b<=circumscribe;b++) { float xcos = cos((b*2*cPI)/circumscribe); float ysin = sin((b*2*cPI)/circumscribe); float xcosr = xcos * radius; float ysinr = ysin * radius; v[0] = vp[0] + vx[0] * xcosr + vy[0] * ysinr; v[1] = vp[1] + vx[1] * xcosr + vy[1] * ysinr; v[2] = vp[2] + vx[2] * xcosr + vy[2] * ysinr; flag=true; i=*(MapLocusEStart(map,v)); if(i) { j=map->EList[i++]; while(j>=0) { ai2 = obj->AtomInfo + cs->IdxToAtm[j]; if((inclH||(!ai2->hydrogen))&& ((!cullByFlag)|| (!(ai2->flags&cAtomFlag_ignore)))) if((j!=a)&&(j!=jj)) { skip_flag=false; if(ai1->vdw==ai2->vdw) { /* handle singularities */ v1 = cs->Coord+3*j; if((v0[0]==v1[0]) && (v0[1]==v1[1]) && (v0[2]==v1[2])) skip_flag=true; } if(ai3->vdw==ai2->vdw) { /* handle singularities */ v1 = cs->Coord+3*j; if((v2[0]==v1[0]) && (v2[1]==v1[1]) && (v2[2]==v1[2])) skip_flag=true; } if(!skip_flag) if(within3f(cs->Coord+3*j,v,ai2->vdw+probe_radius)) { flag=false; break; } } j=map->EList[i++]; } } if(flag&&(dotCnt<stopDot)) { n[0] = vx[0] * xcos + vy[0] * ysin; n[1] = vx[1] * xcos + vy[1] * ysin; n[2] = vx[2] * xcos * vy[2] * ysin; dotCnt++; v+=3; n+=3; I->NDot++; } } } } jj=map2->EList[ii++]; } } } } } MapFree(map2); } MapFree(map); } if((cavity_cull>0)&&(probe_radius>0.75F)&&(!surface_solvent)) { dot_flag=Alloc(int,I->NDot); ErrChkPtr(dot_flag); for(a=0;a<I->NDot;a++) { dot_flag[a]=0; } dot_flag[maxDot]=1; /* this guarantees that we have a valid dot */ map=MapNew(probe_radius_plus,I->Dot,I->NDot,extent); if(map) { MapSetupExpress(map); flag=true; while(flag) { p=dot_flag; v=I->Dot; flag=false; for(a=0;a<I->NDot;a++) { if(!dot_flag[a]) { cnt=0; i=*(MapLocusEStart(map,v)); if(i) { j=map->EList[i++]; while(j>=0) { if(j!=a) { if(within3f(I->Dot+(3*j),v,probe_radius_plus)) { if(dot_flag[j]) { *p=true; flag=true; break; } cnt++; if(cnt>cavity_cull) { *p=true; flag=true; break; } } } j=map->EList[i++]; } } } v+=3; p++; } } MapFree(map); } v = (v0=I->Dot); n = (n0=I->DotNormal); p=dot_flag; c=I->NDot; I->NDot=0; for(a=0;a<c;a++) { if(*(p++)) { *(v0++)=*(v++); *(n0++)=*(n++); *(v0++)=*(v++); *(n0++)=*(n++); *(v0++)=*(v++); *(n0++)=*(n++); I->NDot++; } else { v+=3; n+=3; } } FreeP(dot_flag); } PRINTFD(FB_RepSurface) " GetSolventDots-DEBUG: %d->%d\n",c,I->NDot ENDFD; }