Homec4science
Diffusion lammps (lammm-devel) Browse

lammps/tools/pymol_asphere/src/glsurface.cpp
00aef0fe0014lammm-devel

Last Change
Raw File

glsurface.cpp
View Options

/***************************************************************************
glsurface.cpp
W. Michael Brown
-------------------
Store and manipulate surfaces as OpenGL primitives
begin : Sun Jun 8 2003
copyright : (C) 2003 by W. Michael Brown
email : wmbrown@sandia.gov
***************************************************************************/
#include "glsurface.h"
GLSurface::GLSurface(){
}
GLSurface::~GLSurface(){
}
// Operations on GLline structure
bool operator < (const GLline &one, const GLline &two) {
unsigned minone, mintwo, maxone, maxtwo;
if (one.points[0]<one.points[1]) {
minone=one.points[0];
maxone=one.points[1];
} else {
minone=one.points[1];
maxone=one.points[0];
}
if (two.points[0]<two.points[1]) {
mintwo=two.points[0];
maxtwo=two.points[1];
} else {
mintwo=two.points[1];
maxtwo=two.points[0];
}
if (minone<mintwo)
return true;
else if (minone==mintwo)
if (maxone<maxtwo)
return true;
return false;
}
bool operator == (const GLline &one, const GLline &two) {
if (one.points[0]==two.points[0] && one.points[1]==two.points[1])
return true;
if (one.points[0]==two.points[1] && one.points[1]==two.points[0])
return true;
return false;
}
void GLSurface::clear() {
gllines.clear();
linestrips.clear();
triangles.clear();
vertices.clear();
xmesh.clear();
ymesh.clear();
zmesh.clear();
}
// Reserve vector space for specified number of vertices and triangles
void GLSurface::reserve(unsigned num_v, unsigned num_t) {
vertices.reserve(num_v); triangles.reserve(num_t);
}
// Add a triangle
void GLSurface::addtriangle(Triangle &t) {
triangles.push_back(t);
}
// Add a vertex
void GLSurface::addvertex(Vertex &v) {
vertices.push_back(v);
}
// Add many triangles (via swap)
void GLSurface::addtriangles(vector<Triangle> &t) {
triangles.swap(t);
}
// Add many lines (via swap)
void GLSurface::addlinestrips(vector<LineStrip> &l) {
linestrips.swap(l);
}
// Add a mesh (via swap)
void GLSurface::addxyzmesh(vector<GLline> &x,vector<GLline> &y,
vector<GLline> &z) {
xmesh.swap(x); ymesh.swap(y); zmesh.swap(z);
}
// Preconditions:
// Triangles all loaded with correct vertex ordering
// Vertices all loaded
// Vertices that need to be calculated are marked as invalid
// Postconditions:
// For any normals that are marked as invalid, calculate normals by
// averaging normals for each triangle
void GLSurface::calculatenormals() {
unsigned i;
// Set all invalid normals to zero
for (i=0; i<vertices.size(); i++)
if (vertices[i].valid_normal==false)
vertices[i].normal=vectorPt(0,0,0);
for (i=0; i<triangles.size(); i++) {
vectorPt vec1, vec2, normal;
unsigned id0, id1, id2;
id0=triangles[i].point[0];
id1=triangles[i].point[1];
id2=triangles[i].point[2];
vec1=vertices[id1].cpt-vertices[id0].cpt;
vec2=vertices[id2].cpt-vertices[id0].cpt;
normal[X]=vec1[Z]*vec2[Y]-vec1[Y]*vec2[Z];
normal[Y]=vec1[X]*vec2[Z]-vec1[Z]*vec2[X];
normal[Z]=vec1[Y]*vec2[X]-vec1[X]*vec2[Y];
// Only do vertex averaging if the normal has not been
// explicitly set
if (vertices[id0].valid_normal==false)
vertices[id0].normal+=normal;
if (vertices[id1].valid_normal==false)
vertices[id1].normal+=normal;
if (vertices[id2].valid_normal==false)
vertices[id2].normal+=normal;
}
// Mark all normals as valid
for (i=0; i<vertices.size(); i++)
if (vertices[i].valid_normal==false) {
vertices[i].valid_normal=true;
double no=vertices[i].normal.norm();
if (no==0)
vertices[i].normal=ORIGIN;
else
vertices[i].normal/=no*-1;
}
//error->note[31] << "Calculated normals: " << flush;
return;
}
// Normalize normals.
void GLSurface::normalize() {
unsigned i;
for (i=0; i<vertices.size(); i++) {
vertices[i].normal.normalize();
if (am::not_nan(vertices[i].normal.x)==false) {
cerr << "Bad normal!"; // cerr fix me
vertices[i].normal=vectorPt(0,0,0);
}
}
//error->note[31] << "Normalized normals: " << flush;
}
// Flip normals
void GLSurface::flipnormals() {
unsigned i;
for (i=0; i<vertices.size(); i++)
vertices[i].normal*=-1.0;
//error->note[31] << "Flipped normals: " << flush;
}
// Calculate the Area of the Surface
double GLSurface::surfacearea() {
double sa=0;
for (unsigned i=0; i<triangles.size(); i++)
sa+=triarea(triangles[i]);
return sa;
}
double GLSurface::triarea(Triangle &t) {
double a=vertices[t.point[0]].cpt.dist(vertices[t.point[1]].cpt);
double b=vertices[t.point[1]].cpt.dist(vertices[t.point[2]].cpt);
double c=vertices[t.point[0]].cpt.dist(vertices[t.point[2]].cpt);
double d=a+b+c;
return (sqrt(d*(d-a)*(d-b)*(d-c))/4);
}
// Color surfaces
void GLSurface::color(const colorPt &color) {
unsigned i;
for (i=0; i<vertices.size(); i++)
vertices[i].color=color;
}
void GLSurface::colorbygradient(const colorPt &start, const colorPt &end) {
colorPt step, current;
unsigned i;
step=(end-start)/double(vertices.size());
current=start;
for (i=0; i<vertices.size(); i++) {
vertices[i].color=current;
current+=step;
}
}
void GLSurface::colorbymarker(const colorPt &zero, const colorPt &nonzero) {
unsigned i;
for (i=0; i<vertices.size(); i++) {
if (vertices[i].marker==0)
vertices[i].color=zero;
else
vertices[i].color=nonzero;
}
return;
}
// Color a surface based on interpolation of values on a grid
// void GLSurface::colorbygrid(GridNUM<double> &grid,const colorPt &neg,
// const colorPt &mid, const colorPt &pos, double minv,
// double midv, double maxv) {
// Colors c;
// for (unsigned i=0; i<vertices.size(); i++)
// vertices[i].color=c.gradient(minv,midv,maxv,grid.interp(vertices[i].cpt),
// neg,mid,pos);
// }
void GLSurface::set_transparency(double alpha) {
for (unsigned i=0; i<vertices.size(); i++)
vertices[i].transparency=alpha;
}
// Write out triangles as BEGIN,TRIANGLES,END primitives
void GLSurface::writetris(ofstream &out, const string &objname) {
unsigned i,j,k;
int index;
writepymolheader(out);
out << "BEGIN, TRIANGLES," << endl;
for (i=0; i<triangles.size(); i++) {
// Output coords, colors, and normals of points
for (j=0; j<3; j++) {
index=triangles[i].point[j];
// if (vertices[index].transparency!=1)
out << "ALPHA," << vertices[index].transparency << ",";
out << "COLOR,";
for (k=X; k<=Z; k++)
out << vertices[index].color[k] << ",";
if (vertices[index].valid_normal) {
out << "NORMAL,";
for (k=RED; k<=BLUE; k++)
out << vertices[index].normal[k] << ",";
}
out << "VERTEX,";
for (k=X; k<=Z; k++)
out << vertices[index].cpt[k] << ",";
}
out << endl;
}
out << "END";
writepymoltail(out,objname);
}
// Write out triangles as TRIANGLE primitives
void GLSurface::writetris_surf(ofstream &out, const string &objname) {
unsigned i,j,k;
int index;
writepymolheader(out);
for (i=0; i<triangles.size(); i++) {
// if (vertices[triangles[i].point[0]].transparency!=1)
out << "ALPHA," << vertices[triangles[i].point[0]].transparency << ",";
out << "TRIANGLE,";
// Output coords of points
for (j=0; j<3; j++) {
index=triangles[i].point[j];
for (k=X; k<=Z; k++)
out << vertices[index].cpt[k] << ",";
}
// Output normals of points
for (j=0; j<3; j++) {
index=triangles[i].point[j];
for (k=X; k<=Z; k++)
out << vertices[index].normal[k] << ",";
}
// Output colors of points
for (j=0; j<3; j++) {
index=triangles[i].point[j];
for (k=RED; k<=BLUE; k++)
out << vertices[index].color[k] << ",";
}
out << endl;
}
writepymoltail(out,objname);
}
// Write out triangles as POINTS
void GLSurface::writetris_points(ofstream &out, const string &objname) {
unsigned i,k;
writepymolheader(out);
out << "BEGIN, POINTS," << endl;
for (i=0; i<vertices.size(); i++) {
// Output coords, colors, and normals of points
// if (vertices[i].transparency!=1)
out << "ALPHA," << vertices[i].transparency << ",";
out << "COLOR,";
for (k=X; k<=Z; k++)
out << vertices[i].color[k] << ",";
if (vertices[i].valid_normal) {
out << "NORMAL,";
for (k=RED; k<=BLUE; k++)
out << vertices[i].normal[k] << ",";
}
out << "VERTEX,";
for (k=X; k<=Z; k++)
out << vertices[i].cpt[k] << ",";
out << endl;
}
out << "END\n";
writepymoltail(out,objname);
}
// Write out triangles as POINTS
void GLSurface::write_vspheres(ofstream &out, const string &objname,
double radius) {
unsigned i,k;
writepymolheader(out);
for (i=0; i<vertices.size(); i++) {
// Output coords, colors, and normals of points
if (vertices[i].transparency==0)
continue;
out << "ALPHA," << vertices[i].transparency << ",";
out << "COLOR,";
for (k=X; k<=Z; k++)
out << vertices[i].color[k] << ",";
out << "SPHERE,";
for (k=X; k<=Z; k++)
out << vertices[i].cpt[k] << ",";
out << radius << ",\n";
}
writepymoltail(out,objname);
}
// Write out triangles as POINTS
void GLSurface::writespheres(ofstream &out, const string &objname) {
unsigned i,k;
writepymolheader(out);
for (i=0; i<spheres.size(); i++) {
// Output coords, colors, and normals of points
if (vertices[i].transparency==0)
continue;
out << "ALPHA," << vertices[spheres[i].i].transparency << ",";
out << "COLOR,";
for (k=X; k<=Z; k++)
out << vertices[spheres[i].i].color[k] << ",";
out << "SPHERE,";
for (k=X; k<=Z; k++)
out << vertices[spheres[i].i].cpt[k] << ",";
out << spheres[i].radius << ",\n";
}
writepymoltail(out,objname);
}
// Write out triangles as mesh
// OPTIMIZE TO REMOVE IDENTICAL LINES
void GLSurface::writetris_mesh(ofstream &out, const string &objname) {
unsigned i,j,k;
int index;
writepymolheader(out);
for (i=0; i<triangles.size(); i++) {
out << "BEGIN, LINE_LOOP," << endl;
// Output coords, colors, and normals of points
for (j=0; j<3; j++) {
index=triangles[i].point[j];
// if (vertices[index].transparency!=1)
out << "ALPHA," << vertices[index].transparency << ",";
out << "COLOR,";
for (k=X; k<=Z; k++)
out << vertices[index].color[k] << ",";
if (vertices[index].valid_normal) {
out << "NORMAL,";
for (k=RED; k<=BLUE; k++)
out << vertices[index].normal[k] << ",";
}
out << "VERTEX,";
for (k=X; k<=Z; k++)
out << vertices[index].cpt[k] << ",";
}
out << "END,";
out << endl;
}
writepymoltail(out,objname);
}
// Write out Line Strips (Loops)
void GLSurface::writelinestrips(ofstream &out, const string &objname) {
unsigned i,j,k, size;
int index;
writepymolheader(out);
//out << "LINEWIDTH, 2," << endl;
for (i=0; i<linestrips.size(); i++) {
size=linestrips[i].line.size();
if (size==1)
out << "BEGIN, POINTS," << endl;
else if (size==2)
out << "BEGIN, LINES," << endl;
else if (linestrips[i].loop)
out << "BEGIN, LINE_LOOP," << endl;
else
out << "BEGIN, LINE_STRIP," << endl;
// Output coords, colors, and normals of points
for (j=0; j<linestrips[i].line.size(); j++) {
index=linestrips[i].line[j];
// if (vertices[index].transparency!=1)
out << "ALPHA," << vertices[index].transparency << ",";
out << "COLOR,";
for (k=X; k<=Z; k++)
out << vertices[index].color[k] << ",";
if (vertices[index].valid_normal) {
out << "NORMAL,";
for (k=RED; k<=BLUE; k++)
out << vertices[index].normal[k] << ",";
}
out << "VERTEX,";
for (k=X; k<=Z; k++)
out << vertices[index].cpt[k] << ",";
}
out << "END,";
out << endl;
}
writepymoltail(out,objname);
}
// Add a line
void GLSurface::addline(unsigned v1, unsigned v2) {
GLline line;
line.points[0]=v1;
line.points[1]=v2;
gllines.push_back(line);
}
// Write out GLlines
void GLSurface::writelines(ofstream &out, const string &objname) {
writelines(gllines, out, objname);
}
void GLSurface::writelines(vector<GLline> &l, ofstream &out,
const string &objname) {
unsigned i,j,k;
int index;
writepymolheader(out);
//out << "LINEWIDTH, 2," << endl;
out << "BEGIN, LINES," << endl;
for (i=0; i<l.size(); i++)
for (j=0; j<2; j++) {
index=l[i].points[j];
// if (vertices[index].transparency!=1)
out << "ALPHA," << vertices[index].transparency << ",";
out << "COLOR,";
for (k=X; k<=Z; k++)
out << vertices[index].color[k] << ",";
if (vertices[i].valid_normal) {
out << "NORMAL,";
for (k=RED; k<=BLUE; k++)
out << vertices[index].normal[k] << ",";
}
out << "VERTEX,";
for (k=X; k<=Z; k++)
out << vertices[index].cpt[k] << ",";
out << endl;
}
out << "END,\n";
writepymoltail(out,objname);
}
// Write a set of triangles as an xyzmesh - Requires surface was created on
// a grid of 'resolution' with minimum latice point at 'min'
void GLSurface::writexyzmesh(ofstream &out, const string &objname) {
writelines(xmesh,out,objname+"x");
writelines(ymesh,out,objname+"y");
writelines(zmesh,out,objname+"z");
return;
}
void GLSurface::writepymolheader(ofstream &out) {
out << "from pymol.cgo import *\n"
<< "from pymol import cmd\n"
<< "obj = [\n";
}
void GLSurface::writepymoltail(ofstream &out, const string &objname) {
out << "]\n"
<< "cmd.load_cgo(obj,'" << objname << "')\n\n";
}
// Helper function to set the xyz and the normal for a ellipsoid vertex
void GLSurface::SQE_helper(Vertex &ver, const cPt &rad, const double u,
const double v, const double n,
const double e) {
double sqCun=am::sign(cos(u))*powf(fabs(cos(u)),n);
double sqSun=am::sign(sin(u))*powf(fabs(sin(u)),n);
double sqCve=am::sign(cos(v))*powf(fabs(cos(v)),e);
double sqSve=am::sign(sin(v))*powf(fabs(sin(v)),e);
double sqCu2n=am::sign(cos(u))*powf(fabs(cos(u)),2.0-n);
double sqSu2n=am::sign(sin(u))*powf(fabs(sin(u)),2.0-n);
double sqCv2e=am::sign(cos(v))*powf(fabs(cos(v)),2.0-e);
double sqSv2e=am::sign(sin(v))*powf(fabs(sin(v)),2.0-e);
ver.cpt=rad;
ver.cpt.x*=sqCun*sqCve;
ver.cpt.y*=sqCun*sqSve;
ver.cpt.z*=sqSun;
ver.normal.x=sqCu2n*sqCv2e/rad.x;
ver.normal.y=sqCu2n*sqSv2e/rad.y;
ver.normal.z=sqSu2n/rad.z;
}
void GLSurface::add_super_ellipsoid(const cPt &cen, const cPt &rad,
const Quaternion &q,
const double n, const double e,
const double u1, const double u2,
const double v1, const double v2,
const unsigned u_segs,
const unsigned v_segs,
const colorPt &color, const double alpha){
RotMat rot(q);
Vertex ver;
ver.color=color;
ver.transparency=alpha;
ver.valid_normal=true;
double dU=(u2-u1)/u_segs;
double dV=(v2-v1)/v_segs;
double U=u1;
for (unsigned i=0; i<u_segs; i++) {
double V=v1;
for (unsigned j=0; j<v_segs; j++) {
unsigned offset=vertices.size();
SQE_helper(ver,rad,U,V,n,e);
ver.cpt=rot*ver.cpt+cen;
ver.normal=rot*ver.normal;
addvertex(ver);
SQE_helper(ver,rad,U+dU,V,n,e);
ver.cpt=rot*ver.cpt+cen;
ver.normal=rot*ver.normal;
addvertex(ver);
SQE_helper(ver,rad,U+dU,V+dV,n,e);
ver.cpt=rot*ver.cpt+cen;
ver.normal=rot*ver.normal;
addvertex(ver);
SQE_helper(ver,rad,U,V+dV,n,e);
ver.cpt=rot*ver.cpt+cen;
ver.normal=rot*ver.normal;
addvertex(ver);
Triangle t;
t.point[0]=offset;
t.point[1]=offset+1;
t.point[2]=offset+3;
addtriangle(t);
t.point[0]=offset+1;
t.point[1]=offset+2;
t.point[2]=offset+3;
addtriangle(t);
V+=dV;
}
U+=dU;
}
}
void GLSurface::add_ellipsoid(const cPt &center, const cPt &rad,
const Quaternion &rot, const colorPt &color,
const double alpha, const unsigned res) {
assert(res!=0);
add_super_ellipsoid(center,rad,rot,1.0,1.0,HALFPI*-1.0,HALFPI,PI*-1.0,PI,
res,res,color,alpha);
}
void GLSurface::add_ellipsoid(const cPt &center, const cPt &rad,
const Quaternion &rot, const colorPt &color,
const double alpha) {
add_super_ellipsoid(center,rad,rot,1.0,1.0,HALFPI*-1.0,HALFPI,PI*-1.0,PI,
5,5,color,alpha);
}
c4science · Help