1 #ifndef HAVE_GLUT
main(int argc,char ** argv)2 int main(int argc, char **argv) {
3
4 return 0;
5 }
6 #else
7 /*
8 g++ g++ -I/usr/X11R6/include -L/usr/X11R6/lib -I/Home/romhegg/matthey/packages/glut-3.7/include -L/Home/romhegg/matthey/packages/glut-3.7/lib/glut -Wall xyzviz.cpp -O9 -ffast-math -finline-functions -funroll-loops -DNDEBUG -lglut -lGLU -lGL -lXmu -lXt -lSM -lICE -lXext -lX11 -lXi -lXext -lX11 -lm -o xyzviz
9 g++ -Wall xyzviz.cpp -O9 -ffast-math -finline-functions -funroll-loops -DNDEBUG -I$HOME/include/ -L$HOME/lib/ -L/usr/X11R6/lib -lglut -lGLU -lGL -lm -lXmu -lXi -lX11 -o xyzviz
10 xlC xyzviz.cpp -O2 -I/usr/lpp/OpenGL/include/ -L/usr/lpp/OpenGL/lib/ -lglut -lGLU -lGL -lm -lXmu -lXi -lX11 -o xyzviz
11 */
12
13 #include "DCDTrajectoryReader.h"
14 #include "PPMWriter.h"
15 #include "PGMWriter.h"
16 #include "PDBReader.h"
17 #include "XYZBinReader.h"
18 #include "XYZReader.h"
19 #include "XYZTrajectoryReader.h"
20 #include "mathutilities.h"
21 #include "systemutilities.h"
22 #include "Report.h"
23 #include "Timer.h"
24 #include "openglutilities.h"
25 #include "Matrix3by3.h"
26
27 #include <vector>
28 #include <map>
29 #include <string>
30 #include <iostream>
31 #include <iomanip>
32 #include <fstream>
33 #include <algorithm>
34
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <GL/glut.h>
38
39 using namespace ProtoMol;
40 using namespace ProtoMol::Report;
41 using std::vector;
42 using std::string;
43 using std::map;
44
45 //_____________________________________________________________________ Opaque
46 struct Opaque {
OpaqueOpaque47 Opaque():i(0),x(0.0),y(0.0),z(0.0),o(0.0){}
OpaqueOpaque48 Opaque(int a, Vector3D b, double c):i(a),x(b.x),y(b.y),z(b.z),o(c){}
49 int i;
50 double x,y,z,o;
51 };
52
53 //_____________________________________________________________________ Feedback3Dcolor
54 struct Feedback3Dcolor {
55 GLfloat x;
56 GLfloat y;
57 GLfloat z;
58 GLfloat red;
59 GLfloat green;
60 GLfloat blue;
61 GLfloat alpha;
62 };
63
64 //_____________________________________________________________________ DepthIndex
65 struct DepthIndex {
DepthIndexDepthIndex66 DepthIndex():ptr(NULL),depth(0.0){}
DepthIndexDepthIndex67 DepthIndex(GLfloat *a, GLfloat b):ptr(a),depth(b){}
68
69 GLfloat *ptr;
70 GLfloat depth;
71 };
72
73
74 //_____________________________________________________________________ Enum
75 enum MENU {ANIM_IDLE=-1,
76 PRINT_BUFFER=0,
77 SAVE_FRAMEBUFFER,
78 SAVE_FRAMEBUFFER_EPS,
79 QUIT,
80 TOGGLE_VIEW_CUBE,
81 TOGGLE_PLANE_ROTATION,
82 POLYGON_MODE,
83 TOGGLE_LIGHTING,
84 ANIM_SAVE,
85 ANIM_ON,
86 ANIM_OFF};
87
88 enum POLYMODE {FILL=0,
89 LINE,
90 POINT,
91 LAST}; // Just for modulo
92
93 //_____________________________________________________________________ Forward declaration
94 static bool cmpOpaque (const Opaque& a1, const Opaque& a2);
95
96 static bool nextFrame();
97 static void animationCallback(int value);
98 static void bmovieCallback(int);
99 static void bmovieCallback(long value);
100 static void display(void);
101 static void doAnimation(bool flag);
102 static void frameRateAnimationCallback(int value);
103 static void getFrame(int i);
104 static void keyboardCallback(unsigned char key, int x, int y);
105 static void mainMenuCallback(int value);
106 static void motion(int x, int y);
107 static void mouse(int button, int state, int x, int y);
108 static void parseCmdLine(int argc,
109 char* argv[],
110 string& inputname,
111 string& outputname,
112 double& scale,
113 double& size,
114 int& d,
115 double& from,
116 double& to,
117 Vector3D& rot,
118 int& tracer,
119 int& startFieldOfView,
120 double& fieldOfView,
121 GLfloat& angp,
122 GLfloat& angp2,
123 bool& batch,
124 double& dfr);
125 static void readXYZ(const string& inputname, const Vector3D& rot, int tracer);
126 static void reshape(int w, int h);
127 static void saveFramebufferEPSToFile();
128 static void saveFramebufferToFile();
129 static void startAnimation();
130 static void stopAnimation();
131 static void timerCallback(int value);
132 static void timerCallback(long value);
133 static void timerIdleCallback();
134 static void updateLighting();
135 static void updatePolygonMode();
136 static void usage(char* name);
137 static void verbose();
138 static void visibleParticlesCallback(int value);
139 static double brightness();
140 static Timer timer;
141
142 //_____________________________________________________________________ Statics
143 static bool moving=false;
144 static GLfloat angle = 0; // in degrees
145 static GLfloat angle2 = 0; // in degrees
146 static GLfloat anglePlane = 0; // in degrees
147 static GLfloat anglePlane2 = 0; // in degrees
148 static int startx=0;
149 static int starty=0;
150 static bool movingFieldOfView=false;
151 static int startFieldOfView=0;
152 static double fieldOfView = 80;
153 static double xL=0;
154 static double yL=0;
155 static double zL=0;
156 static double xT=0;
157 static double yT=0;
158 static double zT=0;
159 static double maxL=0;
160 static double scale = 1.0;
161 static double opaque = 0.6;
162 static bool viewCube = true;
163 static bool planeRotation = false;
164 static GLsizei winX = 600;
165 static GLsizei winY = 600;
166 static double radius = 0;
167 static int currentFrame = 0;
168 static POLYMODE polygonMode = FILL;
169 static bool lighting = true;
170 static GLint bufferSize = 0;
171 static bool batch = false;
172 static double dFocusRadius = 0.0;
173 static int subdivisions = 10;
174 static GLfloat lineWidth = 1.0;
175 static GLfloat pointSize = 1.0;
176
177 static bool saveAnimFrame = false;
178 static int timeout = -1;
179 static int frame = 0;
180 static int frameRate = -1;
181 static int frameStep = 1;
182 static int totalFrames = 1;
183 static bool movie = false;
184 static int animation;
185 static int mainMenu;
186 static const int EPSILON_TIME = 5; // [ms]
187 static int frameTimeN = 0;
188 static double from = 0.0;
189 static double to = 0.0;
190 static string outputname = "Noname";
191 static vector<Vector3D> r;
192 static vector<Vector3D> rTmp;
193 static vector<Vector3D> rAll;
194 static vector<string> name;
195 static map<string,Vector3D> colorMap;
196 static map<string,int> color2index;
197 static vector<string> index2color;
198 static map<int,bool> colorVisible;
199 static vector<Opaque> opaqueSort;
200
201 //_____________________________________________________________________ brightness
brightness()202 double brightness(){
203 GLint viewport[4]; // -> x, y, width, height
204 glGetIntegerv(GL_VIEWPORT, viewport);
205 GLint width = viewport[2]-viewport[0];
206 GLint height = viewport[3]-viewport[1];
207
208 long size = width*height*3;
209 GLfloat* pixels = new GLfloat[size];
210 GLint drawBuffer;
211 glGetIntegerv(GL_DRAW_BUFFER, &drawBuffer);
212 glReadBuffer((GLenum)drawBuffer); // read from correct buffer
213 glReadPixels(viewport[0], viewport[1], width, height, GL_RGB, GL_FLOAT, pixels);
214 double count = 0.0;
215 for(long i=0;i<size;i++)
216 count += std::max(std::min((double)pixels[i],1.0),0.0);
217 delete [] pixels;
218 return count/size;
219 }
220
221 //_____________________________________________________________________ cmpOpaque
cmpOpaque(const Opaque & a1,const Opaque & a2)222 bool cmpOpaque (const Opaque& a1, const Opaque& a2){
223 if(a1.o > a2.o)
224 return true;
225 if(a1.o < a2.o)
226 return false;
227 if(a1.z < a2.z)
228 return true;
229 return false;
230 }
231
232 //_____________________________________________________________________ saveFramebufferEPSToFile
saveFramebufferEPSToFile()233 void saveFramebufferEPSToFile(){
234
235 char postFix[20];
236 sprintf(postFix,".%06d.eps",currentFrame);
237 string filename = outputname + string(postFix);
238
239 std::ofstream output(filename.c_str(), std::ios::out);
240 if (!output){
241 report << recoverable << "Attempt to open EPS file \'" << filename
242 << "\' for output failed."<<endr;
243 return;
244 }
245
246 unsigned int count = openglToEPS(output,display);
247
248 output.close();
249 report << plain << "EPS \'" << filename << "\' created containing "<<count<<" OpenGL primitives."<<endr;
250 currentFrame++;
251 }
252
253 //_____________________________________________________________________ updatePolygonMode
updatePolygonMode()254 void updatePolygonMode(){
255 switch (polygonMode) {
256 case POINT:
257 glPolygonMode(GL_FRONT_AND_BACK, GL_POINT);
258 break;
259 case LINE:
260 glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
261 break;
262 case FILL:
263 glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
264 break;
265 case LAST:
266 default:
267 break;
268 }
269 }
270
271 //_____________________________________________________________________ updateLighting
updateLighting()272 void updateLighting(){
273 if (lighting) {
274 glEnable(GL_LIGHTING);
275 } else {
276 glDisable(GL_LIGHTING);
277 }
278 }
279
280
281 //_____________________________________________________________________ saveFramebufferToFile
saveFramebufferToFile()282 void saveFramebufferToFile(){
283
284 char postFix[20];
285 sprintf(postFix,".%06d.ppm",currentFrame);
286 string filename = outputname + string(postFix);
287
288 if(dFocusRadius > 0.0){
289 PGMWriter writer(filename);
290 PGM pgm;
291 openglToPGM(pgm);
292 writer << pgm;
293 }
294 else {
295 PPMWriter writer(filename);
296 PPM ppm;
297 openglToPPM(ppm);
298 writer << ppm;
299 }
300
301 report << plain << "PPM \'" << filename << "\' created, frame "<<frame<<"."<<endr;
302
303 currentFrame++;
304 }
305
306 //_____________________________________________________________________ readXYZ
readXYZ(const string & inputname,const Vector3D & rot,int tracer)307 void readXYZ(const string& inputname, const Vector3D& rot, int tracer){
308
309 XYZ xyz;
310 vector<XYZ> trajectory;
311
312 // PDB
313 if(PDBReader(inputname).tryFormat()){
314 PDBReader reader(inputname);
315 PDB tmp;
316 if(!(reader >> tmp))
317 report << error << "Could not read PDB file \'"<<inputname<<"\'."<<endr;
318 xyz.clear();
319 xyz.coords = tmp.coords;
320 for(unsigned int j=0;j<tmp.atoms.size();++j)
321 xyz.names.push_back(tmp.atoms[j].elementName);
322 trajectory.push_back(xyz);
323 }
324 // XYZ
325 else if(XYZReader(inputname).tryFormat()){
326 XYZReader reader(inputname);
327 if(!(reader >> xyz))
328 report << error << "Could not read XYZ file \'"<<inputname<<"\'."<<endr;
329 trajectory.push_back(xyz);
330 }
331 // Binary XYZ
332 else if(XYZBinReader(inputname).tryFormat()){
333 XYZBinReader reader(inputname);
334 if(!(reader >> xyz))
335 report << error << "Could not read binary XYZ file \'"<<inputname<<"\'."<<endr;
336 trajectory.push_back(xyz);
337 }
338 // DCD
339 else if(DCDTrajectoryReader(inputname).tryFormat()){
340 DCDTrajectoryReader reader(inputname);
341 while((reader >> xyz)){
342 trajectory.push_back(xyz);
343 }
344 //const int n = (argc >1 ? atoi(argv[2]):0);
345 }
346 // Trajectory XYZ
347 else if(XYZTrajectoryReader(inputname).tryFormat()){
348 XYZTrajectoryReader reader(inputname);
349 while((reader >> xyz)){
350 trajectory.push_back(xyz);
351 }
352 }
353 else if(!isAccessible(inputname)){
354 report << error<< "Can not open \'"<<inputname<<"\'."<<endr;
355 }
356 else {
357 report << error<< "Can not figure out format of \'"<<inputname<<"\', skipping."<<endr;
358 }
359
360 movie = (trajectory.size() > 1);
361 totalFrames = trajectory.size();
362
363 vector<Vector3D> colors;
364 if(from < to || dFocusRadius > 0){
365 colors.push_back(Vector3D(0.0,0.0,0.0));
366 }
367 else {
368 colors.push_back(Vector3D(1.0,0.0,0.0));
369 colors.push_back(Vector3D(0.0,1.0,0.0));
370 colors.push_back(Vector3D(0.0,0.0,1.0));
371 colors.push_back(Vector3D(0.0,1.0,1.0));
372 colors.push_back(Vector3D(1.0,1.0,0.0));
373 colors.push_back(Vector3D(1.0,0.0,1.0));
374 colors.push_back(Vector3D(0.0,0.0,0.0));
375 colors.push_back(Vector3D(0.7,0.7,0.7));
376 }
377
378 colorMap.clear();
379 color2index.clear();
380 index2color.clear();
381 colorVisible.clear();
382 r.clear();
383 rAll.clear();
384 Matrix3by3 mat;
385 mat.rotate(rot,Vector3D(0.0,0.0,1.0));
386
387 Vector3D a,b;
388 for(unsigned int i=0;i<trajectory.size();++i){
389 for(unsigned int j=0;j<trajectory[i].size();++j){
390 trajectory[i].coords[j] = mat*trajectory[i].coords[j];
391 Real x = trajectory[i].coords[j].x;
392 Real y = trajectory[i].coords[j].y;
393 Real z = trajectory[i].coords[j].z;
394 string tempstr = trajectory[i].names[j];
395
396 if(i==0&&(j==0||x<a.x)) a.x=x;
397 if(i==0&&(j==0||y<a.y)) a.y=y;
398 if(i==0&&(j==0||z<a.z)) a.z=z;
399 if(i==0&&(j==0||x>b.x)) b.x=x;
400 if(i==0&&(j==0||y>b.y)) b.y=y;
401 if(i==0&&(j==0||z>b.z)) b.z=z;
402
403 if(tracer > 1 && (i + (trajectory[0].size()%tracer)) % tracer == 0)
404 tempstr += "-Tracer";
405
406 name.push_back(tempstr);
407 rAll.push_back(Vector3D(x,y,z));
408 if(colorMap.find(tempstr) == colorMap.end()){
409 int index = colorMap.size() % colors.size();
410 colorMap[tempstr]=colors[index];
411 color2index[tempstr]=index;
412 index2color.push_back(tempstr);
413 colorVisible[index]=true;
414 }
415 }
416 }
417 report << endr;
418
419 xL = b.x-a.x;
420 yL = b.y-a.y;
421 zL = b.z-a.z;
422 xT = a.x;
423 yT = a.y;
424 zT = a.z;
425 maxL = std::max(xL,std::max(yL,zL));
426 for(unsigned int i = 0; i < rAll.size(); i++){
427 rAll[i].x += xT;
428 rAll[i].y += yT;
429 rAll[i].z += zT;
430 }
431
432 r.resize(trajectory[0].size());
433 radius = pow(maxL*maxL*maxL/r.size(),1.0/3.0)/10.0*scale;
434 getFrame(0);
435 }
436
437 //_____________________________________________________________________ parseCmdLine
parseCmdLine(int argc,char * argv[],string & inputname,string & outputname,double & scale,double & size,int & d,double & f,double & t,Vector3D & rot,int & tracer,int & startview,double & fov,GLfloat & angp,GLfloat & angp2,bool & batch,double & dfr)438 void parseCmdLine(int argc,
439 char* argv[],
440 string& inputname,
441 string& outputname,
442 double& scale,
443 double& size,
444 int& d,
445 double& f,
446 double& t,
447 Vector3D& rot,
448 int& tracer,
449 int& startview,
450 double& fov,
451 GLfloat& angp,
452 GLfloat& angp2,
453 bool& batch,
454 double& dfr){
455 if (argc < 2 || (argc >= 2 && (string(argv[1]) =="-h" ||
456 string(argv[1]) =="--help" ))){
457 usage(argv[0]);
458 exit(0);
459 }
460
461 int cur = 1;
462 while (cur<(argc-1) && argv[cur][0]=='-') {
463
464 string str(argv[cur]);
465
466 if (str == "-scale") {
467 scale = atof(argv[++cur]);
468 cur++;
469 continue;
470 }
471
472 if (str == "-plane" && cur<(argc-2)) {
473 f = atof(argv[++cur]);
474 t = atof(argv[++cur]);
475 cur++;
476 continue;
477 }
478
479 if (str == "-view" && cur<(argc-2)) {
480 startview = atoi(argv[++cur]);
481 fov = atof(argv[++cur]);
482 cur++;
483 continue;
484 }
485
486 if (str == "-angplane" && cur<(argc-2)) {
487 angp = atof(argv[++cur]);
488 angp2 = atof(argv[++cur]);
489 cur++;
490 continue;
491 }
492
493 if (str == "-dfradius" && cur<(argc-1)) {
494 dfr = atof(argv[++cur]);
495 cur++;
496 continue;
497 }
498
499 if (str == "-batch") {
500 batch = true;
501 cur++;
502 continue;
503 }
504
505 if (str == "-size") {
506 size = atof(argv[++cur]);
507 cur++;
508 continue;
509 }
510
511 if (str == "-tracer") {
512 tracer = atoi(argv[++cur]);
513 cur++;
514 continue;
515 }
516
517 if (str == "-subdivisions") {
518 d = atoi(argv[++cur]);
519 cur++;
520 continue;
521 }
522
523 if (str == "-z" || str == "-001"){
524 rot = Vector3D(0,0,1);
525 cur++;
526 continue;
527 }
528
529 if (str == "-y" || str == "-010"){
530 rot = Vector3D(0,1,0);
531 cur++;
532 continue;
533 }
534
535 if (str == "-x" || str == "-100"){
536 rot = Vector3D(1,0,0);
537 cur++;
538 continue;
539 }
540
541 if (str == "-xy" || str == "-110"){
542 rot = Vector3D(1,1,0);
543 cur++;
544 continue;
545 }
546
547 if (str == "-xz" || str == "-101"){
548 rot = Vector3D(1,0,1);
549 cur++;
550 continue;
551 }
552
553 if (str == "-yz" || str == "-011"){
554 rot = Vector3D(0,1,1);
555 cur++;
556 continue;
557 }
558
559 if (str == "-xyz" || str == "-111"){
560 rot = Vector3D(1,1,1);
561 cur++;
562 continue;
563 }
564
565 if (str == "-rot") {
566 rot.x = atof(argv[++cur]);
567 rot.y = atof(argv[++cur]);
568 rot.z = atof(argv[++cur]);
569 cur++;
570 continue;
571 }
572
573
574 break;
575
576 }
577
578 inputname = string(argv[cur++]);
579 if (cur < argc){
580 outputname = string(argv[cur]);
581 }
582 }
583
584 //_____________________________________________________________________ usage
usage(char * name)585 void usage(char* name){
586 report << plain << "Usage: " << name << " [-h][--help] [-scale <real>] [-plane <real> <real> [-movie] [-subdivisions <int>] [-size <real>] [-x][-y][-z][-xy][-xz][-yz][-xyz][-rot <vector>][-tracer <int>] inputname [outputname]\n"
587 << " where:\n"
588 << " -h or --help : This help output\n"
589 << " -scale <real> : (optional) sphere scaling factor\n"
590 << " -subdivisions <int> : (optional) number of subdivisions around/anlong the Z axis for sphere\n"
591 << " -size <real> : (optional) point size and line width\n"
592 << " -movie : (optional) for XYZ trajectory input\n"
593 << " -plane <real> <real> : (optional) z-plane projection\n"
594 << " -x : (optional) projection axis\n"
595 << " -y : (optional) projection axis\n"
596 << " -z : (optional) projection axis\n"
597 << " -xy : (optional) projection axis\n"
598 << " -xz : (optional) projection axis\n"
599 << " -yz : (optional) projection axis\n"
600 << " -xyz : (optional) projection axis\n"
601 << " -rot <vector> : (optional) rotate to (0,0,1)\n"
602 << " -tracer <int> : (optional) color every n'th differently\n"
603 << " inputname : filename of input\n"
604 << " outputname : (optional) filename of EPS/PPM output"<<endr;
605 }
606
607 //_____________________________________________________________________ display
display()608 void display() {
609 double fact = dFocusRadius*dFocusRadius/radius/radius;
610 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
611 glPushMatrix();
612 glTranslatef(-xT, -yT, -5 * maxL);
613 glTranslatef(xT, yT, zT);
614 glRotatef(angle2, 1.0, 0.0, 0.0);
615 glRotatef(angle, 0.0, 1.0, 0.0);
616 glTranslatef(-xT, -yT, -zT);
617 glPushMatrix();
618 glTranslatef(xT, yT, zT);
619 if(viewCube){
620 glColor3f(0.0, 0.0, 0.0);
621 glutWireCube(maxL);
622 }
623 glPopMatrix();
624 glEnable(GL_BLEND);
625 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
626 #ifdef NO_OPAQUE_SORT
627 for(unsigned int n = 0; n < r.size(); n++) {
628 if(colorVisible[color2index[name[n]]]){
629 //std::cerr << r[n].x << "," << r[n].y << "," << r[n].z << "," << radius << "\n";
630 if(from >= to && dFocusRadius <= 0.0){
631 glPushMatrix();
632 glTranslatef(r[n].x, r[n].y, r[n].z);
633 glColor4f(colorMap[name[n]].x, colorMap[name[n]].y, colorMap[name[n]].z, opaque);
634 glutSolidSphere(radius, subdivisions, subdivisions);
635 glPopMatrix();
636 }
637 else {
638 double x = r[n].x-xT;
639 double y = r[n].y-yT;
640 double z = r[n].z-zT;
641 double a1=cos(anglePlane/180.0*M_PI);
642 double a2=sin(anglePlane/180.0*M_PI);
643 double b1=cos(anglePlane2/180.0*M_PI);
644 double b2=sin(anglePlane2/180.0*M_PI);
645 double u = a1*x + a2*b1*y -a2*b2*z;
646 double v = -a2*x + a1*b1*y -a1*b2*z;
647 double w = -b2*y + b1*z;
648 x = u+xT;
649 y = v+yT;
650 z = w+zT;
651 if(dFocusRadius > 0.0){
652 double op = 1.0/(1.0+fact*z*z);
653 glPushMatrix();
654 glTranslatef(x, y, z);
655 glColor4f(colorMap[name[n]].x, colorMap[name[n]].y, colorMap[name[n]].z, opaque*op);
656 glutSolidSphere(radius/sqrt(op), subdivisions, subdivisions);
657 glPopMatrix();
658 }
659 else if(z >= from && z <= to){
660 double op = z-from;
661 op = 1.0-fabs(2.0*op/(to-from)-1.0);
662 //std::cerr <<from << "," << to << "," << r[n].z << ","<<op<< "\n";
663 glPushMatrix();
664 glTranslatef(x, y, z);
665 glColor4f(colorMap[name[n]].x, colorMap[name[n]].y, colorMap[name[n]].z, opaque*op);
666 glutSolidSphere(radius, subdivisions, subdivisions);
667 glPopMatrix();
668 }
669 }
670 }
671 }
672 #else
673 if(opaqueSort.size() != r.size())
674 opaqueSort.resize(r.size());
675 if(rTmp.size() != r.size())
676 rTmp.resize(r.size());
677 unsigned int count = 0;
678 for(unsigned int n = 0; n < r.size(); n++) {
679 if(colorVisible[color2index[name[n]]]){
680 if(from >= to && dFocusRadius <= 0.0){
681 opaqueSort[count]= Opaque(n,r[n],1.0);
682 ++count;
683 }
684 else {
685
686 // Edit here
687 double x = r[n].x-xT;
688 double y = r[n].y-yT;
689 double z = r[n].z-zT;
690 double u,v;
691 u = y*cos(anglePlane2/180.0*M_PI) - z*sin(anglePlane2/180.0*M_PI);
692 v = z*cos(anglePlane2/180.0*M_PI) + y*sin(anglePlane2/180.0*M_PI);
693 y = u;
694 z = v;
695 u = x*cos(anglePlane/180.0*M_PI) - z*sin(anglePlane/180.0*M_PI);
696 v = z*cos(anglePlane/180.0*M_PI) + x*sin(anglePlane/180.0*M_PI);
697 x = u;
698 z = v;
699 // double a1=cos(anglePlane/180.0*M_PI);
700 // double a2=sin(anglePlane/180.0*M_PI);
701 // double b1=cos(anglePlane2/180.0*M_PI);
702 // double b2=sin(anglePlane2/180.0*M_PI);
703 // double u = b1*x + b2*y ;
704 // double v = -a1*b2*x + a1*b1*y + a2*z;
705 // double w = a2*b2*x - a2*b1*y + a1*z;
706 // x = u;
707 // y = v;
708 // z = w;
709 double z0 = z;
710 x += xT;
711 y += yT;
712 z += zT;
713 if(dFocusRadius > 0.0) {
714 double op = 1.0/(1.0+fact*z0*z0);
715 if(op > 1e-2){
716 opaqueSort[count]=Opaque(n,Vector3D(x,y,z),op);
717 ++count;
718 }
719 }
720 else if(z0 >= from && z0 <= to){
721 double op = z0-from;
722 op = 1.0-fabs(2.0*op/(to-from)-1.0);
723 opaqueSort[count]=Opaque(n,Vector3D(x,y,z),op);
724 ++count;
725 }
726 }
727 }
728 }
729 std::sort(opaqueSort.begin(),opaqueSort.begin()+count,cmpOpaque);
730 for(unsigned int i = 0; i < count; i++) {
731 int n = opaqueSort[i].i;
732 double op = opaqueSort[i].o*opaque;
733 double rad = radius;
734 if(dFocusRadius > 0.0){
735 rad /= sqrt(opaqueSort[i].o);
736 }
737 glPushMatrix();
738 glTranslatef(opaqueSort[i].x,opaqueSort[i].y,opaqueSort[i].z);
739 glColor4f(colorMap[name[n]].x, colorMap[name[n]].y, colorMap[name[n]].z, op);
740 glutSolidSphere(rad, subdivisions, subdivisions);
741 glPopMatrix();
742 }
743 #endif
744 glDisable(GL_BLEND);
745 glPopMatrix();
746 glutSwapBuffers();
747 }
748
749
750 //_____________________________________________________________________ reshape
reshape(int w,int h)751 void reshape (int w, int h) {
752 glViewport(0, 0, (GLsizei) w, (GLsizei) h);
753 glMatrixMode(GL_MODELVIEW);
754 glLoadIdentity();
755 glMatrixMode(GL_PROJECTION);
756 glLoadIdentity();
757 gluPerspective(fieldOfView, w / static_cast<float>(h), 0.1*maxL, 10*maxL);
758 winX = w;
759 winY = h;
760 }
761
762 //_____________________________________________________________________ mouse
mouse(int button,int state,int x,int y)763 void mouse (int button, int state, int x, int y) {
764 switch (button) {
765 case GLUT_LEFT_BUTTON:
766 if (state == GLUT_DOWN) {
767 movingFieldOfView = true;
768 startFieldOfView = y;
769 }
770 if (state == GLUT_UP) {
771 movingFieldOfView = false;
772 }
773 break;
774 case GLUT_MIDDLE_BUTTON:
775 if (state == GLUT_DOWN) {
776 moving = true;
777 startx = x;
778 starty = y;
779 }
780 if (state == GLUT_UP) {
781 moving = false;
782 }
783 break;
784
785 case GLUT_RIGHT_BUTTON:
786 break;
787 default:
788 break;
789 }
790 }
791
792 //_____________________________________________________________________ motion
motion(int x,int y)793 void motion(int x, int y) {
794 if (moving) {
795 if((from < to || dFocusRadius > 0) && planeRotation){
796 anglePlane = anglePlane + (x - startx);
797 anglePlane2 = anglePlane2 + (y - starty);
798 }
799 else {
800 angle = angle + (x - startx);
801 angle2 = angle2 + (y - starty);
802 }
803 startx = x;
804 starty = y;
805 glutPostRedisplay();
806 }
807 if (movingFieldOfView) {
808 fieldOfView = fieldOfView + (y - startFieldOfView)/10.0;
809 startFieldOfView = y;
810 if(fieldOfView < 10){
811 fieldOfView = 10;
812 }
813 if(fieldOfView > 170){
814 fieldOfView = 170;
815 }
816 //std::cerr << fieldOfView << "\n";
817 reshape(winX,winY);
818 glutPostRedisplay();
819 }
820 }
821
mainMenuCallback(int value)822 void mainMenuCallback(int value){
823 switch (value){
824 case TOGGLE_VIEW_CUBE:
825 viewCube = !viewCube;
826 glutPostRedisplay();
827 break;
828 case TOGGLE_PLANE_ROTATION:
829 planeRotation = !planeRotation;
830 glutSetMenu(mainMenu);
831 glutChangeToMenuEntry(2,const_cast<char*>((string("Toggle Plane Rotation : ")+(planeRotation?string("Yes"):string("No"))).c_str()),TOGGLE_PLANE_ROTATION);
832 break;
833 case QUIT:
834 exit(0);
835 break;
836 case SAVE_FRAMEBUFFER:
837 saveFramebufferToFile();
838 break;
839 case PRINT_BUFFER:
840 openglToPlain(std::cerr,display);
841 break;
842 case SAVE_FRAMEBUFFER_EPS:
843 saveFramebufferEPSToFile();
844
845 break;
846 case POLYGON_MODE:
847 polygonMode = static_cast<POLYMODE>((static_cast<int>(polygonMode) + 1) % static_cast<int>(LAST));
848 updatePolygonMode();
849 glutPostRedisplay();
850 break;
851 case TOGGLE_LIGHTING:
852 lighting = !lighting;
853 updateLighting();
854 glutPostRedisplay();
855 break;
856 default:
857 break;
858 }
859 }
860
861 //_____________________________________________________________________ keyboardCallback
animationCallback(int value)862 void animationCallback(int value) {
863 switch (value) {
864 case ANIM_OFF: // Off
865 stopAnimation();
866 break;
867 case ANIM_SAVE:
868 if (saveAnimFrame){ // no-save frames
869 saveAnimFrame = false;
870 glutSetMenu(animation);
871 glutChangeToMenuEntry(5,"Save frames : No",ANIM_SAVE);
872 }
873 else { // save frames
874 saveAnimFrame = true;
875 glutSetMenu(animation);
876 glutChangeToMenuEntry(5,"Save frames : Yes",ANIM_SAVE);
877 glutPostRedisplay();
878 saveFramebufferToFile();
879 }
880 break;
881 case ANIM_ON:
882 startAnimation();
883 break;
884 default:
885 break;
886 }
887 }
888 //_____________________________________________________________________ stopAnimation
bmovieCallback(long value)889 void bmovieCallback(long value){
890 bmovieCallback(static_cast<int>(value));
891 }
bmovieCallback(int v)892 void bmovieCallback(int v){
893 saveFramebufferToFile();
894 startAnimation();
895 }
896 //_____________________________________________________________________ stopAnimation
stopAnimation()897 void stopAnimation(){
898 timer.stop();
899 if(timer.getTime().getRealTime() > 0 && frameTimeN > 0)
900 report << plain <<"Frames/sec: "<<static_cast<double>(frameTimeN)/timer.getTime().getRealTime()<<"."<<endr;
901 timeout = -1;
902 glutIdleFunc(0);
903 if(batch)
904 exit(0);
905 }
906 //_____________________________________________________________________ startAnimation
startAnimation()907 void startAnimation(){
908 timer.stop();
909 timer.reset();
910 timer.start();
911 frameTimeN = 0;
912 if (frameRate > 0){ // Timeout
913 timeout = (frameRate > EPSILON_TIME) ? frameRate - EPSILON_TIME : 1;
914 glutIdleFunc(0);
915 glutTimerFunc(timeout,&timerCallback,1);
916 }
917 else { // Idle, as fast as possible
918 timeout = 0;
919 glutIdleFunc(timerIdleCallback);
920 }
921 }
922
923 //_____________________________________________________________________ frameRateAnimationCallback
frameRateAnimationCallback(int value)924 void frameRateAnimationCallback(int value){
925 frameRate = value;
926 if(timeout >= 0)
927 startAnimation();
928 }
929
930 //_____________________________________________________________________ visibleParticlesCallback
visibleParticlesCallback(int value)931 void visibleParticlesCallback(int value){
932 colorVisible[value] = !colorVisible[value];
933 display();
934 }
935
936 //_____________________________________________________________________ timerCallback
timerCallback(long value)937 void timerCallback(long value){
938 timerCallback(static_cast<int>(value));
939 }
timerCallback(int value)940 void timerCallback(int value){
941 Timer t;
942 t.start();
943 doAnimation(false);
944 t.stop();
945 int n = static_cast<int>(floor(t.getTime().getRealTime()*1000.0));
946 if (timeout>0){
947 if(timeout-n > 0)
948 glutTimerFunc(timeout-n,timerCallback,value + 1);
949 else
950 glutTimerFunc(1,timerCallback,value + 1);
951 }
952 }
953 //_____________________________________________________________________ timerIdleCallback
timerIdleCallback()954 void timerIdleCallback(){
955 doAnimation(false);
956 }
957
958 //_____________________________________________________________________ doAnimation
doAnimation(bool flag)959 void doAnimation(bool flag){
960 if (timeout >= 0 || flag){
961 if(nextFrame()){
962 display();
963 if (saveAnimFrame)
964 saveFramebufferToFile();
965 }
966 else {
967 if(timeout >= 0){
968 report << plain << "Last frame ("<<frame<<")."<<endr;
969 }
970 stopAnimation();
971 }
972 }
973 else {
974 display();
975 }
976 }
977
978 //_____________________________________________________________________ nextFrame
nextFrame()979 bool nextFrame(){
980 int j = frame;
981 frame += frameStep;
982 int i = frame;
983
984 if(i < 0)
985 i = 0;
986 else if(i >= totalFrames){
987 i = totalFrames -1;
988 }
989
990 if(j < 0)
991 j = 0;
992 else if(j >= totalFrames){
993 j = totalFrames -1;
994 }
995
996 bool update = (i != j);
997 if(update){
998 getFrame(i);
999 frameTimeN++;
1000 }
1001 frame = i;
1002 return update;
1003 }
1004 //_____________________________________________________________________ getFrame
getFrame(int i)1005 void getFrame(int i){
1006 std::copy(&(rAll[i*r.size()]), &(rAll[(i+1)*r.size()]), r.begin());
1007 }
1008
1009 //_____________________________________________________________________ frameStepAnimationCallback
frameStepAnimationCallback(int value)1010 void frameStepAnimationCallback(int value){
1011 frameStep = value;
1012 }
1013
1014 //_____________________________________________________________________ verbose
verbose()1015 void verbose(){
1016 report << plain <<"Particle(s) : " << r.size() << "\n";
1017 report <<"Frame(s) : " << totalFrames << "\n";
1018 report <<"Box : [" << xL << "," << yL << "," << zL << "]\n";
1019 report <<"Types(s) : " << colorMap.size() << "\n";
1020 report <<"Sphere radius : " << radius << "\n";
1021 report <<"Sphere subdivisions : " << subdivisions << "\n";
1022 report <<"Scaling : " << scale << "\n";
1023 report <<"Opaque : " << opaque << "\n";
1024 report <<"Point/line size : " << lineWidth << "\n";
1025 report <<"View cube : " << (viewCube?string("on"):string("off")) << "\n";
1026 report <<"Lighting : " << (lighting?string("on"):string("off")) << "\n";
1027 report <<"Switch fill mode : " << (polygonMode==FILL?string("fill"):(polygonMode==LINE?string("line"):string("point"))) << "\n";
1028 report <<"Window : " << winX<<"x"<<winY<<", angles ("<<angle<<", "<<angle2<<") angles plane ("<<anglePlane<<", "<<anglePlane2<<"), start ("<<startx<<", "<<starty<<"), "
1029 << "view ("<<startFieldOfView<<", "<<fieldOfView<<"), buffer: "<<bufferSize<<".\n";
1030 report <<"Picture number : " << currentFrame << "\n";
1031 report <<"Animation : frame: " << frame <<", frame step: " << frameStep <<", frames/sec:";
1032 if(frameRate > 0)
1033 report <<1000.0/static_cast<double>(frameRate);
1034 else
1035 report <<"max";
1036 if(timer.getActualTime().getRealTime() > 0 && frameTimeN > 0)
1037 std::cerr <<" showing at frames/sec: "<<static_cast<double>(frameTimeN)/timer.getActualTime().getRealTime();
1038 std::cerr <<".\n";
1039 if(from < to) {
1040 report <<"Projection plane : " <<from<<" - "<<to<<endr;
1041 }
1042 else if(dFocusRadius > 0) {
1043 report <<"Projection focus : " <<dFocusRadius<<std::endl;
1044 report <<"Max radius scaling : " <<sqrt(1.0+dFocusRadius*dFocusRadius/radius/radius*maxL*maxL/4.0)<<std::endl;
1045 report <<"Brightness : " <<brightness()<<endr;
1046 }
1047 else {
1048 report <<"3D"<<endr;
1049 }
1050 }
1051 //_____________________________________________________________________ keyboardCallback
keyboardCallback(unsigned char key,int,int)1052 void keyboardCallback(unsigned char key, int, int){
1053 switch (key) {
1054 case 'z':
1055 {
1056 double b = 0.0;
1057 GLfloat a1 = anglePlane;
1058 GLfloat a2 = anglePlane2;
1059 GLfloat b1 = anglePlane;
1060 GLfloat b2 = anglePlane2;
1061
1062 for(int i = -5;i <= 5;i++){
1063 for(int j = -5;j <= 5;j++){
1064 anglePlane = b1+i;
1065 anglePlane2 = b2+j;
1066 display();
1067 double a = brightness();
1068 if(a > b){
1069 a1 = anglePlane;
1070 a2 = anglePlane2;
1071 b = a;
1072 }
1073 }
1074 }
1075 anglePlane = a1;
1076 anglePlane2 = a2;
1077 display();
1078 }
1079 break;
1080 case 'w': // Save framebuffer PPM
1081 saveFramebufferToFile();
1082 break;
1083 case 'p': // Save framebuffer EPS
1084 saveFramebufferEPSToFile();
1085 break;
1086 case 27: // ESC and 'q' will quit.
1087 case 'q':
1088 exit(1);
1089 break;
1090 case '-':
1091 opaque -= 0.05;
1092 if(opaque < 0.0)
1093 opaque = 0.0;
1094 else
1095 glutPostRedisplay();
1096 break;
1097 case '+':
1098 opaque += 0.05;
1099 if(opaque > 1.0)
1100 opaque = 1.0;
1101 else
1102 glutPostRedisplay();
1103 break;
1104 case 'r':
1105 planeRotation = !planeRotation;
1106 glutPostRedisplay();
1107 break;
1108 case ' ':
1109 viewCube = !viewCube;
1110 glutPostRedisplay();
1111 break;
1112 case 'v':
1113 verbose();
1114 break;
1115 case 'a': // Animation statistic
1116 if(movie){
1117 report << plain <<"Frame: " << frame <<", frame step: " << frameStep <<", frames/sec:";
1118 if(frameRate > 0)
1119 report <<1000.0/static_cast<double>(frameRate);
1120 else
1121 report <<"max";
1122 report <<"."<<endr;
1123
1124 }
1125 break;
1126 case 's': // Stop animation
1127 if(movie){
1128 stopAnimation();
1129 }
1130 break;
1131 case 'i': // Idle animation (as fast as possible)
1132 if(movie){
1133 startAnimation();
1134 }
1135 break;
1136
1137 case 'f': // Stop animation and one step forward
1138 if(movie){
1139 stopAnimation();
1140 doAnimation(true);
1141 }
1142 break;
1143
1144 case 'b': // Stop animation and one step back
1145 if(movie){
1146 stopAnimation();
1147 frameStep *= -1;
1148 doAnimation(true);
1149 frameStep *= -1;
1150 }
1151 break;
1152 case 'c': // Stop and reset animation
1153 if(movie){
1154 stopAnimation();
1155 frame = 0; // Go back to frame 0
1156 frameStep = 1; // frameStep = 1
1157 getFrame(frame);
1158 glutPostRedisplay();
1159 }
1160 break;
1161
1162 case 9: // Change sign of frameStep
1163 if(movie){
1164 frameStep *=-1;
1165 }
1166 break;
1167 case 'h':
1168 case '?':
1169 report << plain <<"Summary of keyboard commands:" << std::endl
1170 << " q, ESC : Quit." << std::endl
1171 << " w : Write PPM file." << std::endl
1172 << " p : Write EPS file." << std::endl
1173 << " SPACE : Toggle view cube." << std::endl
1174 << " +, - : Increase/decrease colors." << std::endl
1175 << " h, ? : This message." << std::endl
1176 << " v : Debug information." << std::endl;
1177 if(from < to || dFocusRadius > 0.0)
1178 report << " r : Toggle plane rotation." << std::endl;
1179 if(movie){
1180 report << " s : Stop animation." << std::endl
1181 << " i : Start animation." << std::endl
1182 << " f : Stop animation and one step forward." << std::endl
1183 << " b : Stop animation and one step back." << std::endl
1184 << " a : Animation statistic." << std::endl
1185 << " TAB : Change sign of frameStep." << std::endl
1186 << " c : Set actual frame = 0 and frameStep = 1." << endr;
1187 }
1188
1189 break;
1190 default:
1191 break;
1192 }
1193 }
1194
1195 //_____________________________________________________________________ main
main(int argc,char ** argv)1196 int main (int argc, char** argv) {
1197 glutInit(&argc, argv);
1198
1199
1200 string inputname;
1201 double size = 0.0;
1202 int d = 0;
1203 Vector3D normal(1.0,0,0);
1204 int tracer = 0;
1205 parseCmdLine(argc,argv, inputname, outputname, scale, size, d, from, to, normal,tracer,startFieldOfView,fieldOfView,anglePlane,anglePlane2,batch,dFocusRadius);
1206
1207 if(size > 0){
1208 lineWidth = size;
1209 pointSize = size;
1210 }
1211 if(d > 2){
1212 subdivisions = d;
1213 }
1214
1215 if(batch || (from < to) || (dFocusRadius > 0.0))
1216 viewCube = false;
1217
1218 planeRotation = (from < to) || (dFocusRadius > 0.0);
1219
1220 if(dFocusRadius > 0.0)
1221 opaque = 1.0;
1222
1223 glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
1224 glutInitWindowSize(winX, winY);
1225 glutInitWindowPosition(100, 100);
1226 glutCreateWindow(argv[0]);
1227
1228 GLfloat whiteLight[] = {1.0, 1.0, 1.0, 1.0};
1229 GLfloat lightPosition1[] = {5, 5, 10.0, 0.0};
1230
1231 glClearColor(1.0, 1.0, 1.0, 0.0);
1232 glShadeModel(GL_SMOOTH);
1233
1234 glLightfv(GL_LIGHT0, GL_DIFFUSE, whiteLight);
1235 glLightfv(GL_LIGHT0, GL_POSITION, lightPosition1);
1236 glEnable(GL_LIGHTING);
1237 glEnable(GL_LIGHT0);
1238 glEnable(GL_DEPTH_TEST);
1239 glEnable(GL_COLOR_MATERIAL);
1240 glEnable(GL_CULL_FACE);
1241
1242 glLineWidth(lineWidth);
1243 glPointSize(pointSize);
1244
1245 updateLighting();
1246 updatePolygonMode();
1247
1248 glutMotionFunc(motion);
1249 glutDisplayFunc(display);
1250 glutReshapeFunc(reshape);
1251 glutMouseFunc(mouse);
1252 glutKeyboardFunc(keyboardCallback);
1253
1254 readXYZ(inputname,normal,tracer);
1255
1256 // Animation
1257 if(movie){
1258 int frameStepAnimation = glutCreateMenu(frameStepAnimationCallback);
1259 glutAddMenuEntry("-100", -100);
1260 glutAddMenuEntry(" -50", -50);
1261 glutAddMenuEntry(" -20", -20);
1262 glutAddMenuEntry(" -10", -10);
1263 glutAddMenuEntry(" -5", -5);
1264 glutAddMenuEntry(" -2", -2);
1265 glutAddMenuEntry(" -1", -1);
1266 glutAddMenuEntry(" 1", 1);
1267 glutAddMenuEntry(" 2", 2);
1268 glutAddMenuEntry(" 5", 5);
1269 glutAddMenuEntry(" 10", 10);
1270 glutAddMenuEntry(" 20", 20);
1271 glutAddMenuEntry(" 50", 50);
1272 glutAddMenuEntry(" 100", 100);
1273
1274 int frameRateAnimation = glutCreateMenu(frameRateAnimationCallback);
1275 glutAddMenuEntry("25", 40);
1276 glutAddMenuEntry("20", 50);
1277 glutAddMenuEntry("10", 100);
1278 glutAddMenuEntry(" 5", 200);
1279 glutAddMenuEntry(" 4", 250);
1280 glutAddMenuEntry(" 2", 500);
1281 glutAddMenuEntry(" 1", 1000);
1282 glutAddMenuEntry("as fast as possible",ANIM_IDLE);
1283
1284 animation = glutCreateMenu(animationCallback);
1285 glutAddMenuEntry("Start animation", ANIM_ON);
1286 glutAddMenuEntry("Stop animation", ANIM_OFF);
1287 glutAddSubMenu ("Frames/sec", frameRateAnimation);
1288 glutAddSubMenu ("Framestep", frameStepAnimation);
1289 glutAddMenuEntry("Save frames : No", ANIM_SAVE);
1290 }
1291
1292 int visibleParticles = glutCreateMenu(visibleParticlesCallback);
1293 for(unsigned int i=0;i<index2color.size();i++)
1294 glutAddMenuEntry(const_cast<char*>(index2color[i].c_str()),i);
1295
1296 mainMenu = glutCreateMenu(mainMenuCallback);
1297 glutAddMenuEntry("Toggle view cube",TOGGLE_VIEW_CUBE);
1298 if(from < to || dFocusRadius > 0.0)
1299 glutAddMenuEntry(const_cast<char*>((string("Toggle Plane Rotation : ")+(planeRotation?string("Yes"):string("No"))).c_str()),TOGGLE_PLANE_ROTATION);
1300 glutAddMenuEntry("Toggle lighting", TOGGLE_LIGHTING);
1301 glutAddMenuEntry("Save framebuffer PPM", SAVE_FRAMEBUFFER);
1302 glutAddMenuEntry("Save framebuffer EPS", SAVE_FRAMEBUFFER_EPS);
1303
1304 if(movie){
1305 glutAddSubMenu("Animation",animation);
1306 }
1307 glutAddSubMenu("Toggle particles",visibleParticles);
1308 glutAddMenuEntry("Print framebuffer", PRINT_BUFFER);
1309 glutAddMenuEntry("Switch fill mode (poly, line, point)", POLYGON_MODE);
1310 glutAddMenuEntry("Quit", QUIT);
1311 glutAttachMenu(GLUT_RIGHT_BUTTON);
1312
1313 verbose();
1314
1315 if(batch){
1316 saveAnimFrame = true;
1317 glutSetMenu(animation);
1318 glutChangeToMenuEntry(5,"Save frames : Yes",ANIM_SAVE);
1319 glutPostRedisplay();
1320 glutTimerFunc(1000,&bmovieCallback,1);
1321 }
1322
1323 glutMainLoop();
1324 return 0;
1325 }
1326 #endif
1327