/* 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 */ static const 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} }; static const 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} }; static const 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 *) (void *) 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