Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80670224
glsurface.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Sep 2, 00:14
Size
16 KB
Mime Type
text/x-c
Expires
Wed, Sep 4, 00:14 (2 d)
Engine
blob
Format
Raw Data
Handle
20338478
Attached To
rLAMMPS lammps
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 ¢er, 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 ¢er, 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);
}
Event Timeline
Log In to Comment