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

Sphere.c

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

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

#include"Base.h"
#include"Sphere.h"
#include"Vector.h"
#include"Err.h"

#include"MemoryDebug.h"

/* Twelve vertices of icosahedron on unit sphere */
#define tau 0.8506508084F      /* t=(1+sqrt(5))/2, tau=t/sqrt(1+t^2)  */
#define one 0.5257311121F      /* one=1/sqrt(1+t^2) , unit sphere     */

const static float start_points[13][3] = {
{  tau,  one,    0 },
{ -tau,  one,    0 },
{ -tau, -one,    0 },
{  tau, -one,    0 },
{  one,   0 ,  tau },
{  one,   0 , -tau },
{ -one,   0 , -tau },
{ -one,   0 ,  tau },
{   0 ,  tau,  one },
{   0 , -tau,  one },
{   0 , -tau, -one },
{   0 ,  tau, -one }};

const static int icosahedron[21][3] = {
    { 4, 8, 7 },
    { 4, 7, 9 },
    { 5, 6, 11  },
    { 5, 10, 6 },
    { 0, 4, 3 },
    { 0, 3, 5 },
    { 2, 7, 1 },
    { 2, 1, 6 },
    { 8, 0, 11 },
    { 8, 11, 1 },
    { 9, 10, 3 },
    { 9, 2, 10 },
    { 8, 4, 0 },
    { 11, 0, 5 },
    { 4, 9, 3 },
    { 5, 3, 10 },
    { 7, 8, 1 },
    { 6, 1, 11 },
    { 7, 2, 9 },
    { 6, 10, 2 }
};

const static int mesh[30][2] = {
{ 0 , 3 },
{ 0 , 4 },
{ 0 , 5 },
{ 0 , 8 },
{ 0 , 11 },
{ 1 , 2 },
{ 1 , 6 },
{ 1 , 7 },
{ 1 , 8 },
{ 1 , 11 },
{ 2 , 6 },
{ 2 , 7 },
{ 2 , 9 },
{ 2 , 10 },
{ 3 , 4 },
{ 3 , 5 },
{ 3 , 9 },
{ 3 , 10 },
{ 4 , 7 },
{ 4 , 8 },
{ 4 , 9 },
{ 5 , 6 },
{ 5 , 10 },
{ 5 , 11 },
{ 6 , 10 },
{ 6 , 11 },
{ 7 , 8 },
{ 7 , 9 },
{ 8 , 11 },
{ 9 , 10 }                                                                        
};


#define FAST_SPHERE_INIT


#ifdef FAST_SPHERE_INIT
#include"SphereData.h"

#else

static SphereRec *MakeDotSphere(PyMOLGlobals *G,int level);

#endif

#if 0
static void SphereDump(FILE *f,char *prefix,SphereRec *sp)
{
  int a;
  int c;
  fprintf(f,"static int %s_NTri = %d;\n",prefix,sp->NTri);
  fprintf(f,"static int %s_NStrip = %d;\n",prefix,sp->NStrip);
  fprintf(f,"static int %s_NVertTot = %d;\n",prefix,sp->NVertTot);
  fprintf(f,"static int %s_nDot = %d;\n",prefix,sp->nDot);
  fprintf(f,"static float %s_dot[][3] = {\n",prefix);
  for(a=0;a<sp->nDot;a++) {
    fprintf(f,"{ %15.12fF, %15.12fF, %15.12fF },\n",
            sp->dot[a][0],
            sp->dot[a][1],
            sp->dot[a][2]);
  }
  fprintf(f,"};\n");
  
  fprintf(f,"static float %s_area[] = {\n",prefix);

  c = 0;
  for(a=0;a<sp->nDot;a++) {
    fprintf(f,"%15.12fF,",
            sp->area[a]);
    c = (c+1)%4;
    if(!c) fprintf(f,"\n");
  }
  fprintf(f,"};\n");

  fprintf(f,"static int %s_StripLen[] = {\n",prefix);
  c = 0;
  for(a=0;a<sp->NStrip;a++) {
    fprintf(f,"%6d,",
            sp->StripLen[a]);
    c = (c+1)%10;
    if(!c) fprintf(f,"\n");
  }
  fprintf(f,"};\n");

  fprintf(f,"static int %s_Sequence[] = {\n",prefix);
  c = 0;
  for(a=0;a<sp->NVertTot;a++) {
    fprintf(f,"%6d,",
            sp->Sequence[a]);
    c = (c+1)%10;
    if(!c) fprintf(f,"\n");
  }
  fprintf(f,"};\n");

  fprintf(f,"static int %s_Tri[] = {\n",prefix);
  c = 0;
  for(a=0;a<3*sp->NTri;a++) {
    fprintf(f,"%6d,",
            sp->Tri[a]);
    c = (c+1)%10;
    if(!c) fprintf(f,"\n");
  }
  fprintf(f,"};\n");


}

static void SphereDumpAll(void)
{
  FILE *f;
  f = fopen("SphereData.h","w");
  SphereDump(f,"Sphere0",I->Sphere[0]);
  SphereDump(f,"Sphere1",I->Sphere[1]);
  SphereDump(f,"Sphere2",I->Sphere[2]);
  SphereDump(f,"Sphere3",I->Sphere[3]);
  SphereDump(f,"Sphere4",I->Sphere[4]);
  fclose(f);
}
#endif


void SphereInit(PyMOLGlobals *G)
{
  register CSphere *I = (G->Sphere = Calloc(CSphere,1));

#ifdef FAST_SPHERE_INIT
  I->Array = Alloc(SphereRec,5);

  I->Array[0].area = Sphere0_area;
  I->Array[0].dot = Sphere0_dot;
  I->Array[0].StripLen = Sphere0_StripLen;
  I->Array[0].Sequence = Sphere0_Sequence;
  I->Array[0].NStrip = Sphere0_NStrip;
  I->Array[0].NVertTot = Sphere0_NVertTot;
  I->Array[0].nDot = Sphere0_nDot;
  I->Array[0].Tri = Sphere0_Tri;
  I->Array[0].NTri = Sphere0_NTri;
  I->Array[0].Mesh = (int*)mesh;
  I->Array[0].NMesh = 30;

  I->Array[1].area = Sphere1_area;
  I->Array[1].dot = Sphere1_dot;
  I->Array[1].StripLen = Sphere1_StripLen;
  I->Array[1].Sequence = Sphere1_Sequence;
  I->Array[1].NStrip = Sphere1_NStrip;
  I->Array[1].NVertTot = Sphere1_NVertTot;
  I->Array[1].nDot = Sphere1_nDot;
  I->Array[1].Tri = Sphere1_Tri;
  I->Array[1].NTri = Sphere1_NTri;
  I->Array[1].Mesh = NULL;
  I->Array[1].NMesh = 0;

  I->Array[2].area = Sphere2_area;
  I->Array[2].dot = Sphere2_dot;
  I->Array[2].StripLen = Sphere2_StripLen;
  I->Array[2].Sequence = Sphere2_Sequence;
  I->Array[2].NStrip = Sphere2_NStrip;
  I->Array[2].NVertTot = Sphere2_NVertTot;
  I->Array[2].nDot = Sphere2_nDot;
  I->Array[2].Tri = Sphere2_Tri;
  I->Array[2].NTri = Sphere2_NTri;
  I->Array[2].Mesh = NULL;
  I->Array[2].NMesh = 0;

  I->Array[3].area = Sphere3_area;
  I->Array[3].dot = Sphere3_dot;
  I->Array[3].StripLen = Sphere3_StripLen;
  I->Array[3].Sequence = Sphere3_Sequence;
  I->Array[3].NStrip = Sphere3_NStrip;
  I->Array[3].NVertTot = Sphere3_NVertTot;
  I->Array[3].nDot = Sphere3_nDot;
  I->Array[3].Tri = Sphere3_Tri;
  I->Array[3].NTri = Sphere3_NTri;
  I->Array[3].Mesh = NULL;
  I->Array[3].NMesh = 0;

  I->Array[4].area = Sphere4_area;
  I->Array[4].dot = Sphere4_dot;
  I->Array[4].StripLen = Sphere4_StripLen;
  I->Array[4].Sequence = Sphere4_Sequence;
  I->Array[4].NStrip = Sphere4_NStrip;
  I->Array[4].NVertTot = Sphere4_NVertTot;
  I->Array[4].nDot = Sphere4_nDot;
  I->Array[4].Tri = Sphere4_Tri;
  I->Array[4].NTri = Sphere4_NTri;
  I->Array[4].Mesh = NULL;
  I->Array[4].NMesh = 0;

  I->Sphere[0] = &I->Array[0];
  I->Sphere[1] = &I->Array[1];
  I->Sphere[2] = &I->Array[2];
  I->Sphere[3] = &I->Array[3];
  I->Sphere[4] = &I->Array[4];
#else
  I->Sphere[0] = MakeDotSphere(G,0);
  I->Sphere[1] = MakeDotSphere(G,1);
  I->Sphere[2] = MakeDotSphere(G,2);
  I->Sphere[3] = MakeDotSphere(G,3);
  I->Sphere[4] = MakeDotSphere(G,4);
  /*
  SphereDumpAll();
  */
  
#endif


}

#ifndef FAST_SPHERE_INIT
static void SpherePurge(SphereRec *I)
{
  /* NOTE: S->Mesh is not currently a pointer*/
  mfree(I->dot);
  mfree(I->area);
  mfree(I->StripLen);
  mfree(I->Sequence);
  mfree(I->Tri);
  FreeP(I);
}
#endif


void SphereFree(PyMOLGlobals *G)
{
  register CSphere *I = G->Sphere;

#ifndef FAST_SPHERE_INIT
  SpherePurge(I->Sphere[0]);
  SpherePurge(I->Sphere[1]);
  SpherePurge(I->Sphere[2]);
  SpherePurge(I->Sphere[3]);
  SpherePurge(I->Sphere[4]);
#else
  FreeP(I->Array);
#endif
  FreeP(I);
}

/* private stuff */

#ifndef FAST_SPHERE_INIT

#define MAXDOT 2600
#define MAXTRI 6200

typedef int EdgeCol[MAXDOT]; /* should move these into dynamic storage to save 3MB  mem */
typedef EdgeCol EdgeArray[MAXDOT]; 

typedef int Triangle[3];

typedef struct {

  float *Dot;
  EdgeArray *EdgeRef;
  Triangle *Tri;
  int NDot,NTri;
  
} SphereBuilderRec;


static void MakeVertex(SphereBuilderRec *S,int d1,int d2)
{
  if((*S->EdgeRef)[d1][d2]<0)
       {
            average3f(S->Dot+(3*d1),S->Dot+(3*d2),S->Dot+(3*S->NDot));
            (*S->EdgeRef)[d1][d2]=S->NDot;
            (*S->EdgeRef)[d2][d1]=S->NDot;
            normalize3f(S->Dot+(3*S->NDot));
            S->NDot++;
       }
}

static float SphericalAngle(SphereBuilderRec *S,int d0,int d1,int d2)
{
  Vector3f v1,v2,s1,s2;

  /* map vector onto surface of sphere and measure angle */
  subtract3f(S->Dot+(3*d1),S->Dot+(3*d0),v1);
  subtract3f(S->Dot+(3*d2),S->Dot+(3*d0),v2);
  
  remove_component3f(v1,S->Dot+(3*d0),s1);
  remove_component3f(v2,S->Dot+(3*d0),s2);
  return(get_angle3f(s1,s2));
  
}

static SphereRec *MakeDotSphere(PyMOLGlobals *G,int level)
{
  SphereRec *result;
  int *TriFlag;
  int a,b,c,h,k,l,curTri,n,it;
  float area,sumArea=0.0;
  int nStrip,*q,*s;
  int nVertTot;
  int flag;
  float vt1[3],vt2[3],vt[3];
  SphereBuilderRec SBuild,*S;
  S = &SBuild;

  S->Dot=(float*)mmalloc(sizeof(float)*3*MAXDOT);
  ErrChkPtr(G,S->Dot);
  S->EdgeRef=(EdgeArray*)mmalloc(sizeof(EdgeArray));
  ErrChkPtr(G,S->EdgeRef);
  S->Tri=Alloc(Triangle,MAXTRI);
  ErrChkPtr(G,S->Tri);
  TriFlag=Alloc(int,MAXTRI);
  ErrChkPtr(G,TriFlag);
  
  S->NDot = 12;
  for(a=0;a<S->NDot;a++)
       {
            for(c=0;c<3;c++)
              S->Dot[3*a+c]=start_points[a][c];
            normalize3f(S->Dot+(3*a));
       }

  S->NTri = 20;
  for(a=0;a<S->NTri;a++)
       for(c=0;c<3;c++)
            S->Tri[a][c]=icosahedron[a][c];

  for(a=0;a<MAXDOT;a++)
       for(b=0;b<MAXDOT;b++)
            (*S->EdgeRef)[a][b]=-1;

  if(level>4)
       level=4;
  
  for(c=0;c<level;c++)
       {
            /* create new vertices */
            for(a=0;a<S->NTri;a++)
              {
                   MakeVertex(S,S->Tri[a][0],S->Tri[a][1]);
                   MakeVertex(S,S->Tri[a][1],S->Tri[a][2]);
                   MakeVertex(S,S->Tri[a][0],S->Tri[a][2]);
              }         
            /* create new triangles */
            curTri=S->NTri;
            for(a=0;a<curTri;a++)
              {
                   h=S->Tri[a][0];
                   k=S->Tri[a][1];
                   l=S->Tri[a][2];

                   S->Tri[a][0]=h;
                   S->Tri[a][1]=(*S->EdgeRef)[h][k];
                   S->Tri[a][2]=(*S->EdgeRef)[h][l];
                   
                   S->Tri[S->NTri][0]=k;
                   S->Tri[S->NTri][1]=(*S->EdgeRef)[k][h];
                   S->Tri[S->NTri][2]=(*S->EdgeRef)[k][l];
                   S->NTri++;
                   
                   S->Tri[S->NTri][0]=l;
                   S->Tri[S->NTri][1]=(*S->EdgeRef)[l][h];
                   S->Tri[S->NTri][2]=(*S->EdgeRef)[l][k];
                   S->NTri++;
                   
                   S->Tri[S->NTri][0]=(*S->EdgeRef)[h][k];
                   S->Tri[S->NTri][1]=(*S->EdgeRef)[k][l];
                   S->Tri[S->NTri][2]=(*S->EdgeRef)[l][h];
                   S->NTri++;       
            }
            /*          printf( "MakeDotSphere: Level: %i  S->NTri: %i\n",c, S->NTri); */
       }
  /*  printf(" MakeDotSphere: NDot %i S->NTri %i\n",NDot,S->NTri);*/
  result= Alloc(SphereRec,1);
  ErrChkPtr(G,result);
  result->dot = Alloc(Vector3f,S->NDot);
  ErrChkPtr(G,result->dot);
  result->area = Alloc(float,S->NDot);
  ErrChkPtr(G,result->area);
  result->StripLen = Alloc(int,S->NTri*3);
  ErrChkPtr(G,result->StripLen);
  result->Sequence = Alloc(int,S->NTri*3);
  ErrChkPtr(G,result->Sequence);

  for(a=0;a<S->NDot;a++)
       {
            for(c=0;c<3;c++)
              result->dot[a][c]=*(S->Dot+(3*a+c));
            result->area[a] = 0.0;
       }

  /* fix normals so that v1-v0 x v2-v0 is the correct normal */

  for(a=0;a<S->NTri;a++) {
    subtract3f(result->dot[S->Tri[a][1]],result->dot[S->Tri[a][0]],vt1);
    subtract3f(result->dot[S->Tri[a][2]],result->dot[S->Tri[a][0]],vt2);
    cross_product3f(vt1,vt2,vt);
    if(dot_product3f(vt,result->dot[S->Tri[a][0]])<0.0) { /* if wrong, then interchange */
      it=S->Tri[a][2]; S->Tri[a][2]=S->Tri[a][1]; S->Tri[a][1]=it;
       }
  }    

  for(a=0;a<S->NTri;a++)
       {
            area = (float)(SphericalAngle(S,S->Tri[a][0],S->Tri[a][1],S->Tri[a][2]) +
              SphericalAngle(S,S->Tri[a][1],S->Tri[a][0],S->Tri[a][2]) +
                          SphericalAngle(S,S->Tri[a][2],S->Tri[a][0],S->Tri[a][1]) - cPI);
            /* multiply by r^2 to get area */
            sumArea+=area;
            area/=3.0;
            result->area[S->Tri[a][0]] += area;
            result->area[S->Tri[a][1]] += area;
            result->area[S->Tri[a][2]] += area;
       }

  if(fabs(sumArea - (4*cPI))>0.001) {
    printf(" MakeDotSphere: sumArea: %8.6f which is %8.6f Pi\n",sumArea,sumArea/cPI);
       ErrFatal(G,"MakeDotSphere","Area of sphere does not sum to 4*pi!\n");
  }

  
  for(a=0;a<S->NTri;a++)
       TriFlag[a]=false;

  nStrip=0;
  nVertTot=0;
  s=result->StripLen;
  q=result->Sequence;

  /* tesselate the sphere in a semi-efficient fashion...this could definitely be improved */

  flag = true;
  while(flag) {
       flag=false;
       a=0;
       while(a<S->NTri)
            {
              if(!TriFlag[a]) {
                   flag=true;
                   
                   TriFlag[a]=true;
                   *(q++)=S->Tri[a][0];
                   *(q++)=S->Tri[a][1];
                   *(q++)=S->Tri[a][2];
                   n=3;
                   
                   b=0;
                   while(b<S->NTri)
                        {
                          if(!TriFlag[b]) {
                               if(((S->Tri[b][0]==q[-2])&&(S->Tri[b][1]==q[-1]))||
                                     ((S->Tri[b][0]==q[-1])&&(S->Tri[b][1]==q[-2]))) {
                                    *(q++)=S->Tri[b][2];
                                    TriFlag[b]=true;
                                    b=0;
                                    n++;
                               } else if(((S->Tri[b][0]==q[-2])&&(S->Tri[b][2]==q[-1]))||
                                     ((S->Tri[b][0]==q[-1])&&(S->Tri[b][2]==q[-2]))) {
                                    *(q++)=S->Tri[b][1];
                                    TriFlag[b]=true;
                                    b=0;
                                    n++;
                               } else if(((S->Tri[b][2]==q[-2])&&(S->Tri[b][1]==q[-1]))||
                                     ((S->Tri[b][2]==q[-1])&&(S->Tri[b][1]==q[-2]))) {
                                    *(q++)=S->Tri[b][0];
                                    TriFlag[b]=true;
                                    b=0;
                                    n++;
                               } 
                          }
                          b++;
                        }
                   if(n==3) {
                        q[-3]=S->Tri[a][1];
                        q[-2]=S->Tri[a][2];
                        q[-1]=S->Tri[a][0];
                        
                        b=0;
                        while(b<S->NTri)
                          {
                               if(!TriFlag[b]) {
                                    if(((S->Tri[b][0]==q[-2])&&(S->Tri[b][1]==q[-1]))||
                                          ((S->Tri[b][0]==q[-1])&&(S->Tri[b][1]==q[-2]))) {
                                      *(q++)=S->Tri[b][2];
                                      TriFlag[b]=true;
                                      b=0;
                                      n++;
                                    } else if(((S->Tri[b][0]==q[-2])&&(S->Tri[b][2]==q[-1]))||
                                                       ((S->Tri[b][0]==q[-1])&&(S->Tri[b][2]==q[-2]))) {
                                      *(q++)=S->Tri[b][1];
                                      TriFlag[b]=true;
                                      b=0;
                                      n++;
                                    } else if(((S->Tri[b][2]==q[-2])&&(S->Tri[b][1]==q[-1]))||
                                                       ((S->Tri[b][2]==q[-1])&&(S->Tri[b][1]==q[-2]))) {
                                      *(q++)=S->Tri[b][0];
                                      TriFlag[b]=true;
                                      b=0;
                                      n++;
                                    } 
                               }
                               b++;
                          }
                   }
                   if(n==3) {
                        q[-3]=S->Tri[a][2];
                        q[-2]=S->Tri[a][0];
                        q[-1]=S->Tri[a][1];
                        b=0;
                        while(b<S->NTri)
                          {
                               if(!TriFlag[b]) {
                                    if(((S->Tri[b][0]==q[-2])&&(S->Tri[b][1]==q[-1]))||
                                          ((S->Tri[b][0]==q[-1])&&(S->Tri[b][1]==q[-2]))) {
                                      *(q++)=S->Tri[b][2];
                                      TriFlag[b]=true;
                                      b=0;
                                      n++;
                                    } else if(((S->Tri[b][0]==q[-2])&&(S->Tri[b][2]==q[-1]))||
                                                       ((S->Tri[b][0]==q[-1])&&(S->Tri[b][2]==q[-2]))) {
                                      *(q++)=S->Tri[b][1];
                                      TriFlag[b]=true;
                                      b=0;
                                      n++;
                                    } else if(((S->Tri[b][2]==q[-2])&&(S->Tri[b][1]==q[-1]))||
                                                       ((S->Tri[b][2]==q[-1])&&(S->Tri[b][1]==q[-2]))) {
                                      *(q++)=S->Tri[b][0];
                                      TriFlag[b]=true;
                                      b=0;
                                      n++;
                                    } 
                               }
                               b++;
                          }
                   }
                   *(s++) = n;
                   nVertTot+=n;
                   nStrip++;
              }
              a++;
            }
  }
  mfree(S->Dot);
  mfree(S->EdgeRef);
  mfree(TriFlag);
  result->Tri = (int*)S->Tri;
  result->Tri = Realloc(result->Tri,int,S->NTri*3);
  result->NTri = S->NTri;
  result->StripLen = Realloc(result->StripLen,int,nStrip);
  result->Sequence = Realloc(result->Sequence,int,nVertTot);
  result->dot = Realloc(result->dot,Vector3f,S->NDot);
  result->area = Realloc(result->area,float,S->NDot);
  result->nDot = S->NDot;
  result->NStrip = nStrip;
  result->NVertTot = nVertTot;
  result->Mesh = NULL;
  result->NMesh = 0;
  if (!level) { /* provide mesh for S->Sphere[0] only...rest, to do.*/
    result->Mesh=(int*)mesh;
    result->NMesh=30;
  }
  /*
  q=result->Sequence;
  for(a=0;a<result->NStrip;a++)
       {
            printf("%d:",result->StripLen[a]);  
            for(b=0;b<result->StripLen[a];b++)
              {
                   printf("%d ",*(q++));
              }
            printf("\n");
       }
  */

  return(result);
}

#endif

Generated by  Doxygen 1.6.0   Back to index