/* Basic diffusion and FRAP simulation Jonathan Ward and Francois Nedelec, Copyright EMBL Nov. 2007-2009 To compile on the mac: g++ main.cc -framework GLUT -framework OpenGL -o frap To compile on Windows using Cygwin: g++ main.cc -lopengl32 -lglut32 -o frap To compile on Linux: g++ main.cc -lglut -lopengl -o frap */ #include "params.h" #ifdef __APPLE__ # include <GLUT/glut.h> #else # include "GL/glut.h" #endif #include <list> #include <vector> #include <iostream> #include <math.h> #include "../cmath3d/TriangleMesh.h" #include "../cmath3d_v/TriangleMesh_v.h" //link to the global params extern ParamsClass params; extern ActiveMesh mesh; ///whether the simulation is paused (=1) or not (=0) int pauseSimulation=1; ///how big chunks of steps should be taken between drawing int stepGranularity=1; ///OpenGL window size -> to recalculate from windows coords to scene coords int windowSizeX=0, windowSizeY=0; float uhel_x = 0.0, uhel_y = 0.0; float zoom = 1.0; bool MouseActive = false; int last_x = 0, last_y = 0; bool showAxes=true; bool showFrame=false; bool showBox=true; bool showWholeMesh=true; bool showVertexSurface=false; bool showQuadricSurface=false; bool showTEST=false; ///manage \e count steps of the simulation and then print and draw void step(int count) { while ((params.currTime < params.stopTime) && (count > 0)) { //do some job params.currTime += params.incrTime; --count; } if (params.currTime >= params.stopTime) { std::cout << "END OF THE SIMULATION HAS BEEN REACHED\n"; pauseSimulation=1; } std::cout << "\nnow at time: " << params.currTime << "\n"; //ask GLUT to call the display function glutPostRedisplay(); } ///draws a flat frame around the playground void displayFrame(void) { glColor3f(params.sceneBorderColour.r, params.sceneBorderColour.g, params.sceneBorderColour.b); glLineWidth(1); glBegin(GL_LINE_LOOP); glVertex3f( params.sceneOffset.x , params.sceneOffset.y , params.sceneCentre.z ); glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y , params.sceneCentre.z ); glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y + params.sceneSize.y , params.sceneCentre.z ); glVertex3f( params.sceneOffset.x , params.sceneOffset.y + params.sceneSize.y , params.sceneCentre.z ); glEnd(); } ///draws a frame box around the playground void displayBox(void) { glColor3f(params.sceneBorderColour.r, params.sceneBorderColour.g, params.sceneBorderColour.b); glLineWidth(1); //front square/plane glBegin(GL_LINE_LOOP); glVertex3f( params.sceneOffset.x , params.sceneOffset.y , params.sceneOffset.z ); glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y , params.sceneOffset.z ); glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z ); glVertex3f( params.sceneOffset.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z ); glEnd(); //rear square/plane glBegin(GL_LINE_LOOP); glVertex3f( params.sceneOffset.x , params.sceneOffset.y , params.sceneOffset.z + params.sceneSize.z ); glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y , params.sceneOffset.z + params.sceneSize.z ); glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z + params.sceneSize.z ); glVertex3f( params.sceneOffset.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z + params.sceneSize.z ); glEnd(); //sideways "connectors" glBegin(GL_LINES); glVertex3f( params.sceneOffset.x , params.sceneOffset.y , params.sceneOffset.z ); glVertex3f( params.sceneOffset.x , params.sceneOffset.y , params.sceneOffset.z + params.sceneSize.z ); glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y , params.sceneOffset.z ); glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y , params.sceneOffset.z + params.sceneSize.z ); glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z ); glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z + params.sceneSize.z ); glVertex3f( params.sceneOffset.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z ); glVertex3f( params.sceneOffset.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z + params.sceneSize.z ); glEnd(); } void displayAxes(void) { glLineWidth(2); glBegin(GL_LINES); glColor3f(1.0,0.0,0.0); //red -- x axis glVertex3f(params.sceneCentre.x,params.sceneCentre.y,params.sceneCentre.z); glVertex3f(params.sceneCentre.x+params.sceneSize.x/4.f,params.sceneCentre.y,params.sceneCentre.z); glColor3f(0.0,1.0,0.0); //green -- y axis glVertex3f(params.sceneCentre.x,params.sceneCentre.y,params.sceneCentre.z); glVertex3f(params.sceneCentre.x,params.sceneCentre.y+params.sceneSize.y/4.f,params.sceneCentre.z); glColor3f(0.0,0.0,1.0); //blue -- z axis glVertex3f(params.sceneCentre.x,params.sceneCentre.y,params.sceneCentre.z); glVertex3f(params.sceneCentre.x,params.sceneCentre.y,params.sceneCentre.z+params.sceneSize.z/4.f); glEnd(); glLineWidth(1); } int VertexID=10; int vertexLevels=5; void ActiveMesh::displayVertexAndNeigs(void) { //determine some reasonable number of nearest neighbors std::vector< std::vector<size_t> > neigsLeveled; ulm::getVertexNeighbours(*this,VertexID,vertexLevels,neigsLeveled); //make it flat... std::vector<size_t> neigs; for (unsigned int l=0; l < neigsLeveled.size(); ++l) for (unsigned int i=0; i < neigsLeveled[l].size(); ++i) neigs.push_back(neigsLeveled[l][i]); neigsLeveled.clear(); /* std::vector<size_t> neigs; ulm::getVertexImmediateNeighbours(*this,VertexID,neigs); */ glPointSize(3.0f); glBegin(GL_POINTS); // centre vertex: glColor3f(0.0f,0.0f,1.0f); glVertex3f(Pos[VertexID].x,Pos[VertexID].y,Pos[VertexID].z); // neigs: glColor3f(1.0f,0.0f,0.0f); for (unsigned int i=1; i < neigs.size(); ++i) glVertex3f(Pos[neigs[i]].x,Pos[neigs[i]].y,Pos[neigs[i]].z); glEnd(); } void ActiveMesh::displayTEST(void) { //get surface params float surf_coeff[10]; //CalcQuadricSurface_Taubin(VertexID,surf_coeff); CalcQuadricSurface_sphere(5.0f,params.sceneCentre,surf_coeff); std::cout << "\nsurface: " << surf_coeff[0] << " + " << surf_coeff[1] << "*x + " << surf_coeff[2] << "*y + " << surf_coeff[3] << "*z + " << surf_coeff[4] << "*xy + " << surf_coeff[5] << "*xz + " << surf_coeff[6] << "*yz + " << surf_coeff[7] << "*x*x + " << surf_coeff[8] << "*y*y + " << surf_coeff[9] << "*z*z = 0\n"; glPointSize(2.0f); glBegin(GL_POINTS); glColor3f(0.0f,1.0f,0.0f); //determine triangle vertices const Vector3FC& v1=Pos[VertexID]; std::cout << "v1: (" << v1.x << "," << v1.y << "," << v1.z << ")\n"; Vector3F point=v1; //backup... Vector3F origPoint(point); //adapt the point std::cout << "dist=" << GetClosestPointOnQuadricSurface(point,surf_coeff); //display the point glColor3f(0.0f,1.0f,0.0f); //glVertex3f(point.x,point.y,point.z); float tmp1,tmp2; if (GetPointOnQuadricSurface(v1.x,v1.y,tmp1,tmp2,surf_coeff)) { glVertex3f(v1.x,v1.y,tmp1); glVertex3f(v1.x,v1.y,tmp2); } if (GetPointOnQuadricSurface(v1.x,v1.z,tmp1,tmp2,surf_coeff)) { glVertex3f(v1.x,tmp1,v1.z); glVertex3f(v1.x,tmp2,v1.z); } if (GetPointOnQuadricSurface(v1.y,v1.z,tmp1,tmp2,surf_coeff)) { glVertex3f(tmp1,v1.y,v1.z); glVertex3f(tmp2,v1.y,v1.z); } /* std::cout << " (" << origPoint.x << "," << origPoint.y << "," << origPoint.z << ") -> (" << point.x << "," << point.y << "," << point.z << ")\n"; */ glColor3f(0.0f,1.0f,1.0f); //glVertex3f(origPoint.x,origPoint.y,origPoint.z); glEnd(); } void ActiveMesh::displayQuadricSurface(void) { //zkusit nad dalsimi krivkami, paraboloid? //pripadne udelat jinak dopocitavani na povrch.. //bo to nyni dava docela velke vzdalenosti pro vypocitany povrch //get surface params float surf_coeff[10]; //CalcQuadricSurface_Taubin(VertexID,surf_coeff); CalcQuadricSurface_sphere(5.0f,params.sceneCentre,surf_coeff); std::cout << "\nsurface: " << surf_coeff[0] << " + " << surf_coeff[1] << "*x + " << surf_coeff[2] << "*y + " << surf_coeff[3] << "*z + " << surf_coeff[4] << "*xy + " << surf_coeff[5] << "*xz + " << surf_coeff[6] << "*yz + " << surf_coeff[7] << "*x*x + " << surf_coeff[8] << "*y*y + " << surf_coeff[9] << "*z*z = 0\n"; glPointSize(2.0f); glBegin(GL_POINTS); glColor3f(0.0f,1.0f,0.0f); //do the rendering, i.e., take incident triangles //and sample points from them and project //these points on the obtained surface //list of neigs triangle IDs std::vector<size_t> neigsT; ulm::getVertexIncidentTriangles(*this,VertexID,neigsT); //FOR CYCLE BEGIN //over all triangles //for (unsigned int t=0; t < neigsT.size(); ++t) for (unsigned int t=0; t < 1; ++t) { //determine triangle vertices const Vector3FC& v1=Pos[ID[3*neigsT[t] +0]]; const Vector3FC& v2=Pos[ID[3*neigsT[t] +1]]; const Vector3FC& v3=Pos[ID[3*neigsT[t] +2]]; std::cout << "v1: (" << v1.x << "," << v1.y << "," << v1.z << ")\n"; std::cout << "v2: (" << v2.x << "," << v2.y << "," << v2.z << ")\n"; std::cout << "v3: (" << v3.x << "," << v3.y << "," << v3.z << ")\n"; //sweep (coarsely!) across current triangle for (float c=0.1f; c <= 0.9f; c += 0.1f) for (float b=0.1f; b <= (0.9f-c); b += 0.1f) { float a=1.0f -b -c; Vector3F point=a*v1; point+=b*v2; point+=c*v3; //backup... Vector3F origPoint(point); //adapt the point std::cout << "dist=" << GetClosestPointOnQuadricSurface(point,surf_coeff); //display the point glColor3f(0.0f,1.0f,0.0f); glVertex3f(point.x,point.y,point.z); std::cout << " (" << origPoint.x << "," << origPoint.y << "," << origPoint.z << ") -> (" << point.x << "," << point.y << "," << point.z << ")\n"; glColor3f(0.0f,1.0f,1.0f); glVertex3f(origPoint.x,origPoint.y,origPoint.z); } } //FOR CYCLE END glEnd(); } void ActiveMesh::displayMesh(void) { glColor3f(0.6f,0.6f,0.6f); long unsigned int* id=ID.data(); for (unsigned int i=0; i < ID.size(); i+=3) { //glBegin(GL_TRIANGLES); glBegin(GL_LINE_LOOP); glNormal3f(norm[i/3].x,norm[i/3].y,norm[i/3].z); glVertex3f(Pos[*id].x,Pos[*id].y,Pos[*id].z); ++id; glVertex3f(Pos[*id].x,Pos[*id].y,Pos[*id].z); ++id; glVertex3f(Pos[*id].x,Pos[*id].y,Pos[*id].z); ++id; glEnd(); } } void display(void) { glClear(GL_COLOR_BUFFER_BIT); glLoadIdentity(); glTranslatef(params.sceneCentre.x,params.sceneCentre.y,params.sceneCentre.z); glScalef(zoom,zoom,zoom); Vector3d<float> eye=params.sceneCentre; eye.x+=params.sceneSize.z*sinf(uhel_x)*cosf(uhel_y); eye.y+=params.sceneSize.z*sinf(uhel_y); eye.z+=params.sceneSize.z*cosf(uhel_x)*cosf(uhel_y); gluLookAt(eye.x,eye.y,eye.z, //eye params.sceneCentre.x,params.sceneCentre.y,params.sceneCentre.z, //center //0.0, cosf(-uhel_y), sinf(-uhel_y)); // up 0.0, 1.0, 0.0); // up glEnable(GL_FOG); if (showBox) displayBox(); if (showFrame) displayFrame(); if (showAxes) displayAxes(); glDisable(GL_FOG); if (showWholeMesh) mesh.displayMesh(); if (showVertexSurface) mesh.displayVertexAndNeigs(); if (showQuadricSurface) mesh.displayQuadricSurface(); if (showTEST) mesh.displayTEST(); glFlush(); glutSwapBuffers(); } // called automatically after a certain amount of time void timer(int) { //do simulation and then draw step(stepGranularity); //ask GLUT to call 'timer' again in 'delay' milli-seconds: //in the meantime, you can inspect the drawing :-) if (!pauseSimulation) glutTimerFunc(500, timer, 0); } void GetSceneViewSize(const int w,const int h, float& xVisF, float& xVisT, float& yVisF, float& yVisT) { //what is the desired size of the scene const float xBound=params.sceneSize.x + 2.f*params.sceneOuterBorder.x; const float yBound=params.sceneSize.y + 2.f*params.sceneOuterBorder.y; xVisF = params.sceneOffset.x - params.sceneOuterBorder.x; yVisF = params.sceneOffset.y - params.sceneOuterBorder.y; xVisT = xBound + xVisF; yVisT = yBound + yVisF; if ( (float)w * yBound < (float)h * xBound ) //equals to w/h < x/y { //window is closer to the rectangle than the scene shape yVisF *= ((float)h / float(w)) * (xBound/yBound); yVisT *= ((float)h / float(w)) * (xBound/yBound); } else { xVisF *= ((float)w / float(h)) * (yBound/xBound); xVisT *= ((float)w / float(h)) * (yBound/xBound); } } void printKeysHelp(void) { std::cout << "\nKeys:\n"; std::cout << "ESC or 'q': quit ENTER: enters empty line on a console\n"; std::cout << "' ','n','z': simulation goes/stop at various steps\n"; std::cout << "'r': resets the view\n"; std::cout << "'o','O': zoom in,out\n"; std::cout << "'h': displays this help\n"; std::cout << "'d': displays enabled features,displays status\n"; std::cout << "'m': toggle display of the mesh\n"; std::cout << "'v': toggle display of the vertex supporting region\n"; std::cout << "'V': chooses another vertex\n"; std::cout << "'s': toggle display of the fitted quadratic surface \n"; std::cout << "'t': toggle display of the test/debug function\n"; std::cout << "'a': toggle display of the R,G,B <-> x,y,z axes \n"; std::cout << "'f': toggle display of the frame/plane az z=0\n"; std::cout << "'b': toggle display of the bounding box\n"; std::cout << "'x','X': controls for some testing\n"; std::cout << "'i': not used\n"; } void printDispayedFeauters(void) { if (pauseSimulation) std::cout << "\nsimulation is now paused\n"; else std::cout << "\nsimulation is now in progress\n"; std::cout << "simulation time progress " << params.currTime/params.stopTime*100.f << "%, delta time is " << params.incrTime << " minutes\n"; std::cout << "viewing: x_ang=" << uhel_x << ", y_ang=" << uhel_y << ", zoom=" << zoom << "\n"; std::cout << "Currently displayed features:\n"; } // called when a key is pressed void keyboard(unsigned char key, int mx, int my) { //required for calculating the mouse position within the scene coords float xVisF, yVisF; float xVisT, yVisT; float sx,sy; switch (key) { case 27: //this corresponds to 'escape' case 'q': //quit exit(EXIT_SUCCESS); break; case 13: //Enter key -- just empty lines std::cout << std::endl; break; //control execution of the simulation case ' ': //pause/unpause, refresh every 10 cycles if (params.currTime < params.stopTime) pauseSimulation^=1; if (!pauseSimulation) { stepGranularity=1; //= 1 min timer(0); } break; case 'z': //pause/unpause, refresh every 600 cycles if (params.currTime < params.stopTime) pauseSimulation^=1; if (!pauseSimulation) { stepGranularity=100; //= 10 min timer(0); } break; case 'n': //pause, advance by one frame pauseSimulation=1; step(1); break; //view params case 'r': //set scale 1:1 uhel_x = 0.0; uhel_y = 0.0; zoom = 1.0; glutReshapeWindow(int(params.imgSizeX),int(params.imgSizeY)); glutPostRedisplay(); break; case 'o': zoom -= 0.2f; if (zoom < 0.2f) zoom=0.2f; glutPostRedisplay(); break; case 'O': zoom += 0.2f; glutPostRedisplay(); break; //control hints/help case 'h': //print help printKeysHelp(); break; case 'd': //print what is displayed printDispayedFeauters(); break; //display objects case 'm': showWholeMesh^=true; glutPostRedisplay(); break; case 'v': showVertexSurface^=true; glutPostRedisplay(); break; case 'V': ++VertexID; glutPostRedisplay(); break; case 's': showQuadricSurface^=true; glutPostRedisplay(); break; case 't': showTEST^=true; glutPostRedisplay(); break; //annotation controls case 'a': showAxes^=true; glutPostRedisplay(); break; case 'f': showFrame^=true; glutPostRedisplay(); break; case 'b': showBox^=true; glutPostRedisplay(); break; //some testing controls case 'x': glEnable(GL_FOG); glutPostRedisplay(); break; case 'X': glDisable(GL_FOG); glutPostRedisplay(); break; case 'i': //inspect a cell GetSceneViewSize(windowSizeX,windowSizeY, xVisF,xVisT,yVisF,yVisT); sx=(float)mx /(float)windowSizeX * (xVisT-xVisF) + xVisF; sy=(float)(windowSizeY-my)/(float)windowSizeY * (yVisT-yVisF) + yVisF; xVisF=sx=sy; //just to avoid warning, TODO REMOVE /* for (c=agents.begin(); c != agents.end(); c++) if ((*c)->IsPointInCell(Vector3d<float>(sx,sy,0))) { bool wasSelected=(*c)->isSelected; (*c)->isSelected=true; (*c)->ReportState(); (*c)->isSelected=wasSelected; } */ break; } } // called when a special key is pressed void Skeyboard(int key, int mx, int my) { key=mx=my; //TODO REMOVE mx=key; //TODO REMOVE /* //required for calculating the mouse position within the scene coords float xVisF, yVisF; float xVisT, yVisT; float sx,sy; switch (key) { case GLUT_KEY_LEFT: //rotate cell left by 30deg GetSceneViewSize(windowSizeX,windowSizeY, xVisF,xVisT,yVisF,yVisT); sx=(float)mx /(float)windowSizeX * (xVisT-xVisF) + xVisF; sy=(float)(windowSizeY-my)/(float)windowSizeY * (yVisT-yVisF) + yVisF; for (c=agents.begin(); c != agents.end(); c++) if ((*c)->IsPointInCell(Vector3d<float>(sx,sy,0))) { (*c)->desiredDirection += PI/6.f; if ((*c)->desiredDirection > PI) (*c)->desiredDirection -= 2.f*PI; } break; case GLUT_KEY_RIGHT: //rotate cell left by 30deg GetSceneViewSize(windowSizeX,windowSizeY, xVisF,xVisT,yVisF,yVisT); sx=(float)mx /(float)windowSizeX * (xVisT-xVisF) + xVisF; sy=(float)(windowSizeY-my)/(float)windowSizeY * (yVisT-yVisF) + yVisF; for (c=agents.begin(); c != agents.end(); c++) if ((*c)->IsPointInCell(Vector3d<float>(sx,sy,0))) { (*c)->desiredDirection -= PI/6.f; if ((*c)->desiredDirection < -PI) (*c)->desiredDirection += 2.f*PI; } break; } */ } void mouse(int button, int state, int x, int y) { if (state == GLUT_DOWN) MouseActive = true; else if (state == GLUT_UP) MouseActive = false; last_x=button; //TODO REMOVE last_x = x; last_y = y; } void MouseMotion(int x, int y) { if (MouseActive) { uhel_x -= (float(x - last_x) * 0.01f); uhel_y += (float(y - last_y) * 0.01f); last_x = x; last_y = y; glutPostRedisplay(); } } // called when the window is reshaped, arguments are the new window size void reshape(int w, int h) { glViewport(0, 0, w, h); glMatrixMode(GL_PROJECTION); glLoadIdentity(); //remember the new window size... windowSizeX=w; windowSizeY=h; //what will be the final view port float xVisF, yVisF; float xVisT, yVisT; GetSceneViewSize(w,h,xVisF,xVisT,yVisF,yVisT); //glOrtho(xVisF, xVisT, yVisF, yVisT, params.sceneCentre.z-params.sceneSize.z,params.sceneCentre.z+params.sceneSize.z); glOrtho(xVisF, xVisT, yVisF, yVisT, -1000,1000); glMatrixMode(GL_MODELVIEW); } void initializeGL(void) { int Argc=1; char msg[]="meshSurface GUI"; char *Argv[10]; Argv[0]=msg; glutInit(&Argc, Argv); glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA); //glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH); glutInitWindowSize(800, 800); glutInitWindowPosition(600, 100); glutCreateWindow(Argv[0]); glutKeyboardFunc(keyboard); glutSpecialFunc(Skeyboard); glutDisplayFunc(display); glutReshapeFunc(reshape); glutMouseFunc(mouse); glutMotionFunc(MouseMotion); glClearColor(0.0, 0.0, 0.0, 0.0); //enables alfa-blending for every object/drawing glEnable(GL_BLEND); glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glEnable(GL_POINT_SMOOTH); //glEnable(GL_CULL_FACE); //glEnable(GL_DEPTH_TEST); //enable fog: //glEnable(GL_FOG); glFogi(GL_FOG_MODE,GL_LINEAR); GLfloat fogColor[4]={0.4f*params.sceneBorderColour.r, 0.4f*params.sceneBorderColour.g, 0.4f*params.sceneBorderColour.b, 1.0f}; //has no effect glFogfv(GL_FOG_COLOR,fogColor); //density is not required in linear mode //glFogf(GL_FOG_DENSITY,0.3f); //glHint(GL_FOG_HINT,GL_NICEST); //TODO adjust later glFogf(GL_FOG_START,params.sceneOffset.z); glFogf(GL_FOG_END,params.sceneOffset.z+params.sceneSize.z); } void loopGL(void) { glutMainLoop(); } void closeGL(void) { }