Jiaming Guo / flow3

  #include <stdlib.h>
  #include <string>
  #include <fstream>
  #include <algorithm>
  
  // CUDA include
  #ifdef __CUDACC__
  #include "device_launch_parameters.h"
  #include <cuda.h>
  #include <cuda_runtime_api.h>
  #include "cuda_runtime.h"
  #endif
  
  // OPENGL include
  #include <GL/glut.h>
  #include <GL/freeglut.h>
  
  #include "flow.h"
  
  // STIM include
  #include <stim/visualization/gl_aaboundingbox.h>
  #include <stim/parser/arguments.h>
  #include <stim/visualization/camera.h>
  #include <stim/visualization/colormap.h>
  #include <stim/cuda/cudatools/error.h>
  
  
  //********************parameter setting********************
  // overall parameters
  int vX, vY;
  float dx, dy, dz;										// x, y and z image scaling(units/pixel)
  std::string stackdir = "";								// directory where image stacks will be stored
  stim::arglist args;										// create an instance of arglist
  stim::gl_aaboundingbox<float> bb;						// axis-aligned bounding box object
  stim::camera cam;										// camera object
  unsigned num_edge;										// number of edges in the network
  unsigned num_vertex;									// number of vertex in the network
  std::vector<unsigned> pendant_vertex;					// list of pendant vertex index in GT
  std::vector<std::string> menu_option = { "simulation", "build inlet/outlet", "manufacture" };
  stim::flow<float> flow;									// flow object
  float move_pace;										// camera moving parameter
  float u;												// viscosity
  float rou;												// density
  float max_v;
  float min_v;
  int mods;												// special keyboard input
  std::vector<unsigned char> color;						// velocity color map
  std::vector<int> velocity_bar;							// velocity bar
  
  // hard-coded parameters
  float camera_factor = 1.2f;			// start point of the camera as a function of X and Y size
  float orbit_factor = 0.01f;			// degrees per pixel used to orbit the camera
  float zoom_factor = 10.0f;			// zooming factor
  float border_factor = 20.0f;		// border
  float radii_factor = 1.0f;			// radii changing factor
  GLint subdivision = 20;				// slices and stacks

  float default_radius = 5.0f;			// default radii of network vertex

  float delta = 0.01f;				// small discrepancy
  float eps = 20.0f;					// epsilon threshold
  float max_pressure = 0.0f;			// maximum pressure that the channel can bear

  float height_threshold = 100.0f;	// connection height constraint
  float fragment_ratio = 0.0f;		// fragment ratio

  
  // glut event parameters
  int mouse_x;						// window x-coordinate
  int mouse_y;						// window y-coordinate
  bool LTbutton = false;				// true means down while false means up				
  
  // simulation parameters
  bool simulation = false;			// flag indicates simulation mode
  bool color_bound = false;			// flag indicates velocity color map bound
  bool to_select_pressure = false;	// flag indicates having selected a vertex to modify pressure
  unsigned pressure_index;			// the index of vertex that is clicked
  
  // build inlet/outlet parameters
  bool build_inlet_outlet = false;	// flag indicates building inlets and outlets
  bool modified_bridge = false;		// flag indicates having modified inlet/outlet connection

  bool hilbert_curve = false;			// flag indicates enabling hilbert curves constructions
  bool change_fragment = false;		// flag indicates changing fragment for square wave connections

  
  // manufacture parameters
  bool manufacture = false;			// flag indicates manufacture mode
  
  
  //********************helper function*********************
  // get the network basic information
  inline void get_background() {
  
  	pendant_vertex = flow.get_boundary_vertex();
  	num_edge = flow.edges();
  	num_vertex = flow.vertices();
  
  	// set the initial radii
  	flow.init(num_edge, num_vertex);			// initialize flow object
  	for (unsigned i = 0; i < num_edge; i++)

  		flow.set_r(i, default_radius);

  }
  
  // convert from window coordinates to world coordinates
  inline void window_to_world(GLdouble &x, GLdouble &y, GLdouble &z) {
  
  	GLint    viewport[4];
  	GLdouble modelview[16];
  	GLdouble projection[16];
  	GLdouble winX, winY;
  	GLfloat  winZ;
  
  	glGetIntegerv(GL_VIEWPORT, viewport);
  	glGetDoublev(GL_MODELVIEW_MATRIX, modelview);
  	glGetDoublev(GL_PROJECTION_MATRIX, projection);
  
  	winX = (GLdouble)mouse_x;
  	winY = viewport[3] - (GLdouble)mouse_y;
  	glReadPixels((GLint)winX, (GLint)winY, (GLsizei)1, (GLsizei)1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ);
  	gluUnProject(winX, winY, winZ, modelview, projection, viewport, &x, &y, &z);
  }
  
  
  //********************simulation function**********************
  // initialize flow object
  void flow_initialize() {
  

  	flow.set = true;

  	stim::vec3<float> center = bb.center();
  
  	for (unsigned i = 0; i < pendant_vertex.size(); i++) {
  		if (flow.get_vertex(pendant_vertex[i])[0] <= center[0])
  			flow.P[pendant_vertex[i]] = max_pressure - i * delta;	// should set minor discrepancy
  		else
  			flow.P[pendant_vertex[i]] = (i + 1) * delta;			// algorithm treat 0 as no initial pressure
  	}
  }
  
  // find the stable flow state
  void flow_stable_state() {
  	
  	flow.solve_flow(u);
  	flow.get_color_map(max_v, min_v, color, pendant_vertex);
  	color_bound = true;
  
  	velocity_bar.resize(num_edge);
  	for (unsigned i = 0; i < num_edge; i++)
  		velocity_bar[i] = i;
  	std::sort(velocity_bar.begin(), velocity_bar.end(), [&](int x, int y) {return abs(flow.v[x]) < abs(flow.v[y]); });
  }
  
  
  //********************glut function********************
  // dynamically set menu
  // @param num: number of current menu options
  // @param range: range of option to be set from menu_option list
  void glut_set_menu(int num, int range) {
  
  	// remove last time menu options
  	for (int i = 1; i < num + 1; i++)
  		glutRemoveMenuItem(1);
  
  	// set new menu options
  	std::string menu_name;
  	for (int i = 1; i < range + 1; i++) {
  		menu_name = menu_option[i - 1];
  		glutAddMenuEntry(menu_name.c_str(), i);
  	}
  }
  
  // set up the squash transform to whole screen
  void glut_projection() {
  
  	glMatrixMode(GL_PROJECTION);					// load the projection matrix for editing
  	glLoadIdentity();								// start with the identity matrix
  	vX = glutGet(GLUT_WINDOW_WIDTH);				// use the whole screen for rendering
  	vY = glutGet(GLUT_WINDOW_HEIGHT);
  	glViewport(0, 0, vX, vY);						// specify a viewport for the entire window
  	float aspect = (float)vX / (float)vY;			// calculate the aspect ratio
  	gluPerspective(60, aspect, 0.1, 1000000);		// set up a perspective projection
  }
  
  // translate camera to origin
  void glut_modelview() {
  
  	glMatrixMode(GL_MODELVIEW);						// load the modelview matrix for editing
  	glLoadIdentity();								// start with the identity matrix
  	stim::vec3<float> eye = cam.getPosition();		// get the camera position (eye point)
  	stim::vec3<float> focus = cam.getLookAt();		// get the camera focal point
  	stim::vec3<float> up = cam.getUp();				// get the camera "up" orientation
  
  	gluLookAt(eye[0], eye[1], eye[2], focus[0], focus[1], focus[2], up[0], up[1], up[2]);	// set up the OpenGL camera
  }
  
  // glut render function
  void glut_render() {
  
  	glEnable(GL_DEPTH_TEST);
  	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

  	glClearColor(1.0f, 1.0f, 1.0f, 1.0f);

  	glut_projection();
  	glut_modelview();
  
  	if (!simulation && !build_inlet_outlet || manufacture) {
  		flow.glCylinder0();
  	}
  	else {
  		flow.bounding_box();
  		flow.glSolidSphere(max_pressure, subdivision);
  		flow.mark_vertex();
  		flow.glSolidCone(subdivision);
  		flow.glSolidCylinder(subdivision, color);
  		flow.glSolidCuboid();
  	}
  
  	if (build_inlet_outlet) {
  		flow.line_bridge();
  	}
  	
  	if (manufacture) {
  		flow.glSolidCuboid();
  		flow.tube_bridge(subdivision);
  	}
  	
  
  	// render bars
  	// bring up a pressure bar on left
  	if (to_select_pressure) {
  		
  		glMatrixMode(GL_PROJECTION);									// set up the 2d viewport for mode text printing
  		glPushMatrix();
  		glLoadIdentity();
  		vX = glutGet(GLUT_WINDOW_WIDTH);								// get the current window width
  		vY = glutGet(GLUT_WINDOW_HEIGHT);								// get the current window height
  		glViewport(0, 0, vX, vY);										// locate to left bottom corner
  		gluOrtho2D(0, vX, 0, vY);										// define othogonal aspect
  
  		glMatrixMode(GL_MODELVIEW);
  		glPushMatrix();
  		glLoadIdentity();
  
  		glLineWidth(border_factor);
  		glBegin(GL_LINES);
  		glColor3f(0.0, 0.0, 1.0);										// blue to red
  		glVertex2f(border_factor, border_factor);
  		glColor3f(1.0, 0.0, 0.0);
  		glVertex2f(border_factor, (vY - 2.0f * border_factor));
  		glEnd();
  		glFlush();
  
  		// pressure bar text

  		glColor3f(0.0f, 0.0f, 0.0f);

  		glRasterPos2f(0.0f, vY - border_factor);
  		std::stringstream ss_p;
  		ss_p << "Pressure Bar";
  		glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_p.str().c_str()));
  
  		// pressure range text
  		float step = vY - 3.0f * border_factor;
  		step /= 10;
  		for (unsigned i = 0; i < 11; i++) {
  			glRasterPos2f((border_factor * 1.5f), (border_factor + i * step));
  			std::stringstream ss_n;
  			ss_n << (float)i * max_pressure / 10;

  			glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_n.str().c_str()));

  		}
  		glPopMatrix();
  		glMatrixMode(GL_PROJECTION);
  		glPopMatrix();
  	}
  
  	// bring up a velocity bar on left
  	if (simulation && !to_select_pressure) {
  		
  		glMatrixMode(GL_PROJECTION);									// set up the 2d viewport for mode text printing
  		glPushMatrix();
  		glLoadIdentity();
  		vX = glutGet(GLUT_WINDOW_WIDTH);								// get the current window width
  		vY = glutGet(GLUT_WINDOW_HEIGHT);								// get the current window height
  		glViewport(0, 0, vX, vY);										// locate to left bottom corner
  		gluOrtho2D(0, vX, 0, vY);										// define othogonal aspect
  
  		glMatrixMode(GL_MODELVIEW);
  		glPushMatrix();
  		glLoadIdentity();
  
  		float step = (vY - 3 * border_factor);
  		step /= num_edge;
  		for (unsigned i = 0; i < num_edge; i++) {
  			glLineWidth(border_factor);
  			glBegin(GL_LINES);
  			glColor3f((float)color[velocity_bar[i] * 3 + 0] / 255, (float)color[velocity_bar[i] * 3 + 1] / 255, (float)color[velocity_bar[i] * 3 + 2] / 255);
  			glVertex2f(border_factor, border_factor + i * step);
  			glVertex2f(border_factor, border_factor + (i + 1) * step);
  			glEnd();
  		}
  		glFlush();
  
  		// pressure bar text

  		glColor3f(0.0f, 0.0f, 0.0f);

  		glRasterPos2f(0.0f, vY - border_factor);
  		std::stringstream ss_p;
  		ss_p << "Velocity range";
  		glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_p.str().c_str()));
  
  		// pressure range text
  		step = vY - 3 * border_factor;
  		step /= 10;
  		for (unsigned i = 0; i < 11; i++) {
  			glRasterPos2f(border_factor * 1.5f, border_factor + i * step);
  			std::stringstream ss_n;
  			ss_n << min_v + i * (max_v - min_v) / 10;

  			glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_n.str().c_str()));

  		}
  		glPopMatrix();
  		glMatrixMode(GL_PROJECTION);
  		glPopMatrix();
  	}

  
  	// bring up a ratio bar on the left
  	if (change_fragment) {
  		
  		glMatrixMode(GL_PROJECTION);									// set up the 2d viewport for mode text printing
  		glPushMatrix();
  		glLoadIdentity();
  		vX = glutGet(GLUT_WINDOW_WIDTH);								// get the current window width
  		vY = glutGet(GLUT_WINDOW_HEIGHT);								// get the current window height
  		glViewport(0, 0, vX, vY);										// locate to left bottom corner
  		gluOrtho2D(0, vX, 0, vY);										// define othogonal aspect
  
  		glMatrixMode(GL_MODELVIEW);
  		glPushMatrix();
  		glLoadIdentity();
  
  		glLineWidth(border_factor);
  		glBegin(GL_LINES);
  		glColor3f(0.0, 0.0, 1.0);										// blue to red
  		glVertex2f(border_factor, border_factor);
  		glColor3f(1.0, 0.0, 0.0);
  		glVertex2f(border_factor, (vY - 2.0f * border_factor));
  		glEnd();
  		glFlush();
  
  		// pressure bar text
  		glColor3f(0.0f, 0.0f, 0.0f);
  		glRasterPos2f(0.0f, vY - border_factor);
  		std::stringstream ss_p;
  		ss_p << "Ratio bar";
  		glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_p.str().c_str()));
  
  		// pressure range text
  		float step = vY - 3.0f * border_factor;
  		step /= 10;
  		for (unsigned i = 0; i < 11; i++) {
  			glRasterPos2f((border_factor * 1.5f), (border_factor + i * step));
  			std::stringstream ss_n;
  			ss_n << (float)i * 1.0f / 10;
  			glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_n.str().c_str()));
  		}
  		glPopMatrix();
  		glMatrixMode(GL_PROJECTION);
  		glPopMatrix();
  	}
  

  	glutSwapBuffers();
  }
  
  // register glut menu options
  void glut_menu(int value) {
  
  	int num = glutGet(GLUT_MENU_NUM_ITEMS);
  	if (value == 1) {
  		simulation = true;
  		build_inlet_outlet = false;
  		manufacture = false;
  		modified_bridge = false;

  		if (!flow.set)
  			flow_initialize();

  		flow_stable_state();					// main function of solving the linear system
  		flow.print_flow();
  		
  		glut_set_menu(num, 2);
  	}
  
  	if (value == 2) {
  		simulation = false;
  		build_inlet_outlet = true;
  		manufacture = false;
  		if (!modified_bridge) {
  			flow.set_main_feeder();

  			flow.build_synthetic_connection(u, default_radius);

  		}
  
  		glut_set_menu(num, 3);
  	}
  
  	if (value == 3) {
  		simulation = false;
  		build_inlet_outlet = false;
  		manufacture = true;

  		flow.check_synthetic_connection(u, default_radius);

  	}
  
  	glutPostRedisplay();
  }
  
  // defines camera motion based on mouse dragging
  void glut_motion(int x, int y) {
  
  	mods = glutGetModifiers();
  	if (LTbutton && mods == 0) {
  
  		float theta = orbit_factor * (mouse_x - x);		// determine the number of degrees along the x-axis to rotate
  		float phi = orbit_factor * (y - mouse_y);		// number of degrees along the y-axis to rotate
  
  		cam.OrbitFocus(theta, phi);						// rotate the camera around the focal point
  	}
  	mouse_x = x;										// update the mouse position
  	mouse_y = y;
  
  	glutPostRedisplay();								// re-draw the visualization
  }
  
  // get click window coordinates
  void glut_mouse(int button, int state, int x, int y) {
  
  	mouse_x = x;
  	mouse_y = y;
  	if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN)
  		LTbutton = true;
  	else if (button == GLUT_LEFT_BUTTON && state == GLUT_UP)
  		LTbutton = false;
  
  	if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && simulation && !to_select_pressure) {
  
  		GLdouble posX, posY, posZ;
  		window_to_world(posX, posY, posZ);			// get the world coordinates
  
  		bool flag = flow.epsilon_vertex((float)posX, (float)posY, (float)posZ, eps, pressure_index);
  		if (flag) {
  			std::vector<unsigned>::iterator it = std::find(pendant_vertex.begin(), pendant_vertex.end(), pressure_index);
  			if (it != pendant_vertex.end()) 		// if it is dangle vertex
  				to_select_pressure = true;
  		}
  	}
  	else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && simulation && to_select_pressure) {
  		if (y > 2 * border_factor || y < vY - border_factor) {		// within the pressure bar range
  			to_select_pressure = false;
  			float tmp_pressure = (float)(vY - y - border_factor) / ((float)vY - border_factor) * max_pressure;
  			flow.set_pressure(pressure_index, tmp_pressure);
  
  			flow_stable_state();									// main function of solving the linear system
  			flow.print_flow();
  		}
  	}

  	else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && modified_bridge && change_fragment) {
  		if (y > 2 * border_factor || y < vY - border_factor) {		// within the ratio bar range
  			fragment_ratio = (float)(vY - y - border_factor) / ((float)vY - border_factor) * 1.0f;
  			flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold, fragment_ratio, default_radius);
  			change_fragment = false;
  		}
  	}

  }
  
  // define camera move based on mouse wheel move
  void glut_wheel(int wheel, int direction, int x, int y) {
  	
  	mouse_x = x;
  	mouse_y = y;
  
  	GLdouble posX, posY, posZ;
  	window_to_world(posX, posY, posZ);			// get the world coordinates
  
  	if (!to_select_pressure) {
  		bool flag = flow.epsilon_vertex((float)posX, (float)posY, (float)posZ, eps, pressure_index);
  		if (flag && simulation) {
  			float tmp_r;
  			if (direction > 0) {				// increase radii
  				tmp_r = flow.get_radius(pressure_index);
  				tmp_r += radii_factor;
  			}
  			else {
  				tmp_r = flow.get_radius(pressure_index);
  				tmp_r -= radii_factor;
  				if (tmp_r <= 0)

  					tmp_r = default_radius;

  			}
  			flow.set_radius(pressure_index, tmp_r);
  			flow_stable_state();
  			flow.print_flow();
  		}
  		else {
  			if (direction > 0)								// if it is button 3(up), move closer
  				move_pace = zoom_factor;
  			else											// if it is button 4(down), leave farther
  				move_pace = -zoom_factor;
  
  			cam.Push(move_pace);
  		}
  	}
  	
  	glutPostRedisplay();
  }
  
  // define keyboard inputs
  void glut_keyboard(unsigned char key, int x, int y) {
  
  	// register different keyboard operation
  	switch (key) {
  
  		// zooming
  	case 'w':						// if keyboard 'w' is pressed, then move closer
  		move_pace = zoom_factor;
  		cam.Push(move_pace);
  		break;
  	case 's':						// if keyboard 's' is pressed, then leave farther
  		move_pace = -zoom_factor;
  		cam.Push(move_pace);
  		break;
  
  		// output image stack
  	case 'm':
  		if (manufacture) {
  #ifdef __CUDACC__
  			flow.make_image_stack(dx, dy, dz, stackdir);
  #else
  			std::cout << "You need to have a gpu to make image stack, sorry." << std::endl;
  #endif
  		}
  		else if (build_inlet_outlet && !modified_bridge) {

  			modified_bridge = true;

  
  			if (hilbert_curve)
  				flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold);
  			else
  				change_fragment = true;

  		}
  		break;
  	}

  	glutPostRedisplay();
  }
  
  // glut initialization
  void glut_initialize() {
  
  	int myargc = 1;
  	char* myargv[1];
  	myargv[0] = strdup("generate_network_network");
  
  	glutInit(&myargc, myargv);
  	glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA);
  	glutInitWindowPosition(100, 100);							// set the initial window position
  	glutInitWindowSize(1000, 1000);
  	glutCreateWindow("3D flow simulation");
  
  	glutDisplayFunc(glut_render);
  	glutMouseFunc(glut_mouse);
  	glutMotionFunc(glut_motion);
  	glutMouseWheelFunc(glut_wheel);
  	glutKeyboardFunc(glut_keyboard);
  
  	glutCreateMenu(glut_menu);					// create a menu object 
  	glut_set_menu(0, 1);
  	glutAttachMenu(GLUT_RIGHT_BUTTON);			// register right mouse to open menu option
  
  	stim::vec3<float> c = bb.center();			// get the center of the network bounding box
  	// place the camera along the z-axis at a distance determined by the network size along x and y
  	cam.setPosition(c + stim::vec<float>(0, 0, camera_factor * std::max(bb.size()[0], bb.size()[1])));
  	cam.LookAt(c[0], c[1], c[2]);
  }
  
  int main(int argc, char* argv[]) {
  	
  	// add arguments
  	args.add("help", "prints the help");
  	args.add("network", "load network from .obj or .swc file");
  	args.add("maxpress", "maximum allowed pressure in g / units / s^2, default 2 is for blood when units = um", "2", "real value > 0");
  	args.add("viscosity", "set the viscosity of the fluid (in g / units / s), default .00001 is for blood when units = um", ".00001", "real value > 0");
  	args.add("rou", "set the desity of the fluid (in g / units^3), default 1.06*10^-12 is for blood when units = um", ".00000000000106", "real value > 0");

  	args.add("hilbert", "enable hilbert curves connections", "0", "value 1 for enablement");

  	args.add("stackres", "spacing between pixel samples in each dimension(in units/pixel)", "1 1 1", "real value > 0");
  	args.add("stackdir", "set the directory of the output image stack", "", "any existing directory (ex. /home/name/network)");
  	
  	args.parse(argc, argv);								// parse the command line
  
  	// load network
  	if (args["network"].is_set()) {						// load network from user 
  		std::vector<std::string> tmp = stim::parser::split(args["network"].as_string(), '.');
  		if ("obj" == tmp[1])
  			flow.load_obj(args["network"].as_string());
  		else if ("swc" == tmp[1])
  			flow.load_swc(args["network"].as_string());
  		else {
  			std::cout << "Invalid file type" << std::endl;
  			std::exit(1);
  		}
  	}
  	get_background();
  
  	// blood pressure in capillaries range from 15 - 35 torr
  	// 1 torr = 133.3 Pa
  	max_pressure = args["maxpress"].as_float();
  
  	// normal blood viscosity range from 4 - 15 mPa·s(cP)
  	// 1 Pa·s = 1 g / mm / s
  	u = args["viscosity"].as_float();			// g / units / s
  
  	// normally the blood density in capillaries: 1060 kg/m^3 = 1.06*10^-12 g/um^3
  	rou = args["rou"].as_float();
  

  	// check whether to enable hilbert curves or not
  	hilbert_curve = args["hilbert"].as_int();
  

  	// get the vexel and image stack size
  	dx = args["stackres"].as_float(0);
  	dy = args["stackres"].as_float(1);
  	dz = args["stackres"].as_float(2);
  
  	// get the save directory of image stack
  	if (args["stackdir"].is_set())
  		stackdir = args["stackdir"].as_string();
  
  	// glut main loop
  	bb = flow.boundingbox();
  	glut_initialize();
  	glutMainLoop();
  }