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>

  #include <stim/grids/image_stack.h>

  
  
  //********************parameter setting********************
  // overall parameters

  std::string units;										// units used in this program

  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", "adjustment" };

  stim::flow<float> flow;									// flow object

  stim::flow<float> backup;								// flow backup

  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

  float length = 40.0f;									// cuboid length

  float scale = 1.0f;										// scale factor

  bool image_stack = false;								// flag indicates an image stack been loaded
  stim::image_stack<unsigned char, float> S;				// image stack
  float binary_threshold = 128;							// threshold for binary transformation
  float in = 0.0f;										// total input volume flow rate
  float out = 0.0f;
  float Rt = 0.0f;										// total resistance
  float Qn = 0.0f;										// required input total volume flow rate

  GLint dlist;											// simulation display list
  bool undo = false;										// delete display list

  
  // 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

  int picked_x;						// picked window x-coordinate
  int picked_y;						// picked window y-coordinate

  bool LTbutton = false;				// true means down while false means up		

  
  // simulation parameters

  bool render_direction = false;		// flag indicates rendering flow direction for one edge

  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

  bool mark_index = true;				// flag indicates marking the index near the vertex

  bool glyph_mode = false;			// flag indicates rendering glyph for flow velocity field
  bool frame_mode = false;			// flag indicates rendering filament framing structrues

  bool subdivided = false;			// flag indicates subdivision status

  unsigned pressure_index;			// the index of vertex that is clicked

  unsigned direction_index = UINT_MAX;// the index of edge that is pointed at

  std::vector<stim::vec3<float> > back_vertex;	// vertex back up for marking indices

  
  // 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

  bool picked_connection = false;		// flag indicates picked one connection
  bool render_new_connection = false;	// flag indicates rendering new line connection in trasparency
  bool redisplay = false;				// flag indicates redisplay rendering
  bool connection_done = false;		// flag indicates finishing connections

  bool render_flow_rate = false;		// flag indicates rendering total volume flow rate

  unsigned connection_index = UINT_MAX;// the index of connection that is picked

  unsigned port_index = 0;			// inlet (0) or outlet (1)
  stim::vec3<float> tmp_v1, tmp_v2;	// temp vertex
  int coef;							// computational coefficient factor

  
  // manufacture parameters
  bool manufacture = false;			// flag indicates manufacture mode
  
  
  //********************helper function*********************
  // get the network basic information
  inline void get_background() {
  

  	pendant_vertex = flow.get_pendant_vertex();

  	num_edge = flow.edges();
  	num_vertex = flow.vertices();
  
  	// set the initial radii
  	flow.init(num_edge, num_vertex);			// initialize flow object

  	
  	// if no radius information laoded
  	if (!flow.get_radius(0, 0))
  		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);
  }
  

  // convert current image stack into a binary mask
  #ifdef __CUDACC__
  template <typename T, typename F>

  __global__ void binary_transform(size_t N, T* ptr, F threshold) {

  	size_t ix = blockDim.x * blockIdx.x + threadIdx.x;

  	if (ix >= N) return;					// avoid seg-fault
  
  	if (ptr[ix] >= threshold)				// binary transformation

  		ptr[ix] = 0;

  	else
  		ptr[ix] = 255;

  }
  #endif
  

  
  //********************simulation function**********************
  // initialize flow object
  void flow_initialize() {
  

  	flow.set = true;

  	stim::vec3<float> center = bb.center();

  	flow.P.clear();
  	flow.P.resize(num_vertex, 0);									// clear up initialized pressure

  
  	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]); });
  }
  

  // adjustment on input volume flow rate and corresponding flow simulation
  void adjustment() {
  
  	system("CLS");							// clear up console box
  	std::cout << "Please enter the input total volume flow rate: " << std::endl;
  	std::cin >> Qn;
  
  	flow.adjust(in, out, Rt, Qn, u);
  }
  

  
  //********************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) {

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

  		flow.glCylinder0(scale, undo);

  	}

  	else {	
  		flow.bounding_box();		// bounding box
  		if (num_vertex > 100) {						// if the network is big enough (say 100), use display list
  			if (undo) {								// undo rendering list
  				undo = false;
  				glDeleteLists(dlist, 1);
  			}							
  			if (!glIsList(dlist)) {
  				dlist = glGenLists(1);
  				glNewList(dlist, GL_COMPILE);
  				// render network
  				if (!glyph_mode) {
  					flow.glSolidSphere(max_pressure, subdivision, scale);
  					if (mark_index)

  						flow.mark_vertex(back_vertex, scale);

  					//flow.glSolidCone(subdivision);
  					flow.glSolidCylinder(direction_index, color, subdivision, scale);
  				}
  				// render glyphs
  				else
  					flow.glyph(color, subdivision, scale, frame_mode);
  				
  				glEndList();
  			}
  			glCallList(dlist);

  		}

  		else {										// small network
  			// render network
  			if (!glyph_mode) {
  				flow.glSolidSphere(max_pressure, subdivision, scale);
  				if (mark_index) {

  					flow.mark_vertex(back_vertex, scale);

  					//flow.mark_edge();
  				}
  				//flow.glSolidCone(subdivision);
  				flow.glSolidCylinder(direction_index, color, subdivision, scale);
  			}
  			// render glyphs
  			else
  				flow.glyph(color, subdivision, scale, frame_mode);
  		}
  		flow.glSolidCuboid(subdivision, manufacture, length);		// render bus source

  		if (render_direction && !glyph_mode)						// render the flow direction of the vertex pointed

  			flow.glSolidCone(direction_index, subdivision, scale);

  	}
  

  	if (build_inlet_outlet)

  		flow.line_bridge(redisplay);

  	if (manufacture) {

  		flow.glSolidCuboid(subdivision, manufacture, length);
  		flow.tube_bridge(redisplay, subdivision, scale);

  	}

  
  	if (picked_connection && render_new_connection) {
  		glEnable(GL_BLEND);
  		glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
  		glColor4f(0.0f, 0.0f, 0.0f, 0.4f);
  		glBegin(GL_LINE_STRIP);
  		if (!port_index) {
  			glVertex3f(flow.inlet[connection_index].V[1][0], flow.inlet[connection_index].V[1][1], flow.inlet[connection_index].V[1][2]);
  			glVertex3f(tmp_v1[0], tmp_v1[1], tmp_v1[2]);
  			glVertex3f(tmp_v2[0], tmp_v2[1], tmp_v2[2]);
  			glVertex3f(flow.inlet[connection_index].V[2][0], flow.inlet[connection_index].V[2][1], flow.inlet[connection_index].V[2][2]);
  		}
  		else {
  			glVertex3f(flow.outlet[connection_index].V[1][0], flow.outlet[connection_index].V[1][1], flow.outlet[connection_index].V[1][2]);
  			glVertex3f(tmp_v1[0], tmp_v1[1], tmp_v1[2]);
  			glVertex3f(tmp_v2[0], tmp_v2[1], tmp_v2[2]);
  			glVertex3f(flow.outlet[connection_index].V[2][0], flow.outlet[connection_index].V[2][1], flow.outlet[connection_index].V[2][2]);
  		}
  		glEnd();
  		glFlush();
  		glDisable(GL_BLEND);
  	}

  	// 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 || build_inlet_outlet) && !to_select_pressure && !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();
  
  		float step = (vY - 3 * border_factor);

  		step /= BREWER_CTRL_PTS - 1;
  		for (unsigned i = 0; i < BREWER_CTRL_PTS - 1; i++) {

  			glLineWidth(border_factor);
  			glBegin(GL_LINES);

  			glColor3f(BREWERCP[i * 4 + 0], BREWERCP[i * 4 + 1], BREWERCP[i * 4 + 2]);

  			glVertex2f(border_factor, border_factor + i * step);

  			glColor3f(BREWERCP[(i + 1) * 4 + 0], BREWERCP[(i + 1) * 4 + 1], BREWERCP[(i + 1) * 4 + 2]);

  			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();
  

  		// ratio 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()));
  

  		// ratio 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();
  	}
  

  	if (build_inlet_outlet)

  		if (render_flow_rate)
  			flow.display_flow_rate(in, out);

  	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;

  		render_flow_rate = false;

  		manufacture = false;
  		modified_bridge = false;

  		change_fragment = false;

  		connection_done = false;

  		// first time

  		if (!flow.set) {						// only first time simulation called "simulation", ^_^

  			get_background();					// get the graph information

  			back_vertex = flow.back_vertex();	// vertex back up for marking indices

  			flow_initialize();					// initialize flow condition

  			menu_option[0] = "resimulation";
  		}
  
  		// simulation / resimulation

  		flow_stable_state();					// main function of solving the linear system
  		flow.print_flow();
  		

  		if (!glyph_mode)

  			glut_set_menu(num, 2);

  	}
  
  	if (value == 2) {
  		simulation = false;
  		build_inlet_outlet = true;
  		manufacture = false;

  		if (!modified_bridge && !connection_done) {

  			flow.set_main_feeder();

  			flow.build_synthetic_connection(u, default_radius);

  			flow.check_direct_connection();	// check whether direct connections intersect each other
  			connection_done = true;
  		}
  		else if (modified_bridge) {
  			modified_bridge = false;
  			redisplay = true;
  			flow.clear_synthetic_connection();

  		}
  

  		glut_set_menu(num, 4);

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

  		glyph_mode = false;			// manufacuture mode doesn't need flow direction

  		redisplay = true;

  	}
  

  	if (value == 4) {
  		simulation = true;
  		build_inlet_outlet = false;

  		render_flow_rate = false;

  		manufacture = false;
  		
  		adjustment();					// adjust network flow accordingly
  
  		glut_set_menu(num, 1);
  	}
  

  	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
  }
  

  // defines passive mouse motion function
  void glut_passive_motion(int x, int y) {
  

  	mods = glutGetModifiers();

  
  	// check whether the mouse point near to an edge
  	GLdouble posX, posY, posZ;
  	window_to_world(posX, posY, posZ);			// get the world coordinates
  

  	if (simulation || build_inlet_outlet && !mods) {

  		bool flag = flow.epsilon_edge((float)posX, (float)posY, (float)posZ, eps, direction_index);

  		if (flag && !glyph_mode)

  			render_direction = true;

  		else if (!flag && !glyph_mode) {

  			if (render_direction)				// if the direction is displaying currently, do a short delay

  				Sleep(300);

  			render_direction = false;
  			direction_index = -1;
  		}

  		undo = true;

  	}
  

  	if (mods == GLUT_ACTIVE_SHIFT && picked_connection) {
  		render_new_connection = true;

  		size_t i;

  		if (!port_index) {
  			tmp_v1 = stim::vec3<float>(flow.inlet[connection_index].V[1][0], flow.inlet[connection_index].V[1][1] + (float)(picked_y - y), flow.inlet[connection_index].V[1][2]);
  			tmp_v2 = stim::vec3<float>(flow.inlet[connection_index].V[2][0], flow.inlet[connection_index].V[2][1] + (float)(picked_y - y), flow.inlet[connection_index].V[2][2]);
  			i = flow.inlet[connection_index].V.size();
  			if (coef * tmp_v1[1] < coef * flow.inlet[connection_index].V[i - 1][1]) {
  				tmp_v1[1] = flow.inlet[connection_index].V[i - 1][1];
  				tmp_v2[1] = flow.inlet[connection_index].V[i - 1][1];
  			}
  		}
  		else {
  			tmp_v1 = stim::vec3<float>(flow.outlet[connection_index].V[1][0], flow.outlet[connection_index].V[1][1] + (float)(picked_y - y), flow.outlet[connection_index].V[1][2]);
  			tmp_v2 = stim::vec3<float>(flow.outlet[connection_index].V[2][0], flow.outlet[connection_index].V[2][1] + (float)(picked_y - y), flow.outlet[connection_index].V[2][2]);
  			i = flow.outlet[connection_index].V.size();
  			if (coef * tmp_v1[1] < coef * flow.outlet[connection_index].V[i - 1][1]) {
  				tmp_v1[1] = flow.outlet[connection_index].V[i - 1][1];
  				tmp_v2[1] = flow.outlet[connection_index].V[i - 1][1];
  			}
  		}	
  	}
  	else if (mods == GLUT_ACTIVE_CTRL && picked_connection) {
  		render_new_connection = true;
  		if (!port_index) {
  			tmp_v1 = stim::vec3<float>(flow.inlet[connection_index].V[0][0] + (float)(x - picked_x), flow.inlet[connection_index].V[0][1], flow.inlet[connection_index].V[0][2]);
  			tmp_v2 = stim::vec3<float>(flow.inlet[connection_index].V[1][0] + (float)(x - picked_x), flow.inlet[connection_index].V[1][1], flow.inlet[connection_index].V[1][2]);
  			if (tmp_v1[0] < flow.main_feeder[port_index][0] - length / 2) {
  				tmp_v1[0] = flow.main_feeder[port_index][0] - length / 2;
  				tmp_v2[0] = flow.main_feeder[port_index][0] - length / 2;
  			}
  			else if (tmp_v1[0] > flow.main_feeder[port_index][0] + length / 2) {
  				tmp_v1[0] = flow.main_feeder[port_index][0] + length / 2;
  				tmp_v2[0] = flow.main_feeder[port_index][0] + length / 2;
  			}
  		}
  		else {
  			tmp_v1 = stim::vec3<float>(flow.outlet[connection_index].V[0][0] + (float)(x - picked_x), flow.outlet[connection_index].V[0][1], flow.outlet[connection_index].V[0][2]);
  			tmp_v2 = stim::vec3<float>(flow.outlet[connection_index].V[1][0] + (float)(x - picked_x), flow.outlet[connection_index].V[1][1], flow.outlet[connection_index].V[1][2]);
  			if (tmp_v1[0] > flow.main_feeder[port_index][0] + length / 2) {
  				tmp_v1[0] = flow.main_feeder[port_index][0] + length / 2;
  				tmp_v2[0] = flow.main_feeder[port_index][0] + length / 2;
  			}
  			else if (tmp_v1[0] < flow.main_feeder[port_index][0] - length / 2) {
  				tmp_v1[0] = flow.main_feeder[port_index][0] - length / 2;
  				tmp_v2[0] = flow.main_feeder[port_index][0] - length / 2;
  			}
  		}
  	}
  	else
  		render_new_connection = false;
  
  	mouse_x = x;
  	mouse_y = y;
  
  	glutPostRedisplay();							// re-draw the visualization

  }
  

  // get click window coordinates
  void glut_mouse(int button, int state, int x, int y) {
  

  	mods = glutGetModifiers();						// get special keyboard input
  

  	mouse_x = x;
  	mouse_y = y;

  	if (!mods) {
  		picked_connection = false;
  		render_new_connection = false;
  	}

  	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 && !mods && 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, scale, 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 && !mods && 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 - 3.0f * 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 && !mods && 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 - 3.0f * border_factor) * 1.0f;
  		else if (y < 2 * border_factor)
  			fragment_ratio = 1.0f;
  		else if (y > vY - border_factor)
  			fragment_ratio = 0.0f;
  
  		change_fragment = false;

  		render_flow_rate = true;

  		flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold, in, out, 2, fragment_ratio, default_radius);

  	}

  	// move connections along y-axis
  	else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && mods == GLUT_ACTIVE_SHIFT && !modified_bridge && !picked_connection) {	
  		GLdouble posX, posY, posZ;
  		window_to_world(posX, posY, posZ);			// get the world coordinates
  
  		bool flag = flow.epsilon_edge((float)posX, (float)posY, (float)posZ, eps, connection_index, port_index);
  		if (flag) {
  			picked_connection = true;
  			picked_x = x;
  			picked_y = y;
  			if (!port_index)
  				if (flow.inlet[connection_index].V[2][1] > flow.main_feeder[port_index][1])
  					coef = 1;
  				else
  					coef = -1;
  			else
  				if (flow.outlet[connection_index].V[2][1] > flow.main_feeder[port_index][1])
  					coef = 1;
  				else
  					coef = -1;
  		}
  		else
  			picked_connection = false;
  	}
  	else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && mods == GLUT_ACTIVE_SHIFT && !modified_bridge && render_new_connection) {
  		float l = 0.0f;
  		std::vector<typename stim::vec3<float> > V;

  		size_t i;

  		if (!port_index) {
  			i = flow.inlet[connection_index].V.size();
  			if (tmp_v2[1] != flow.inlet[connection_index].V[i - 1][1]) {
  				V.resize(4);
  				V[0] = flow.inlet[connection_index].V[0];
  				V[1] = tmp_v1;
  				V[2] = tmp_v2;
  				V[3] = flow.inlet[connection_index].V[i - 1];
  				std::swap(flow.inlet[connection_index].V, V);
  			}
  			else {
  				V.resize(3);
  				V[0] = flow.inlet[connection_index].V[0];
  				V[1] = tmp_v1;
  				V[2] = tmp_v2;
  				std::swap(flow.inlet[connection_index].V, V);
  			}
  			// calculate new length
  			for (unsigned i = 0; i < flow.inlet[connection_index].V.size() - 1; i++) {
  				l += (flow.inlet[connection_index].V[i + 1] - flow.inlet[connection_index].V[i]).len();
  			}
  			flow.inlet[connection_index].l = l;
  		}
  		else {
  			i = flow.outlet[connection_index].V.size();
  			if (tmp_v2[1] != flow.outlet[connection_index].V[i - 1][1]) {
  				V.resize(4);
  				V[0] = flow.outlet[connection_index].V[0];
  				V[1] = tmp_v1;
  				V[2] = tmp_v2;
  				V[3] = flow.outlet[connection_index].V[i - 1];
  				std::swap(flow.outlet[connection_index].V, V);
  			}
  			else {
  				V.resize(3);
  				V[0] = flow.outlet[connection_index].V[0];
  				V[1] = tmp_v1;
  				V[2] = tmp_v2;
  				std::swap(flow.outlet[connection_index].V, V);
  			}
  			// calculate new length
  			for (unsigned i = 0; i < flow.outlet[connection_index].V.size() - 1; i++) {
  				l += (flow.outlet[connection_index].V[i + 1] - flow.outlet[connection_index].V[i]).len();
  			}
  			flow.outlet[connection_index].l = l;
  		}
  
  		redisplay = true;
  		render_new_connection = false;
  		picked_connection = false;
  
  		flow.check_direct_connection();
  		flow.backup();								// back up direct synthetic connections
  	}
  	// move connections along x-axis
  	else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && mods == GLUT_ACTIVE_CTRL && !modified_bridge && !picked_connection) {
  		GLdouble posX, posY, posZ;
  		window_to_world(posX, posY, posZ);			// get the world coordinates
  
  		bool flag = flow.epsilon_edge((float)posX, (float)posY, (float)posZ, eps, connection_index, port_index);
  		if (flag) {
  			picked_connection = true;
  			picked_x = x;
  			picked_y = y;
  			if (!port_index)
  				coef = 1;
  			else
  				coef = -1;
  		}
  		else
  			picked_connection = false;
  	}
  	else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && mods == GLUT_ACTIVE_CTRL && !modified_bridge && render_new_connection) {
  		float l = 0.0f;
  		if (!port_index) {
  			flow.inlet[connection_index].V[0] = tmp_v1;
  			flow.inlet[connection_index].V[1] = tmp_v2;
  			// calculate new length
  			for (unsigned i = 0; i < flow.inlet[connection_index].V.size() - 1; i++) {
  				l += (flow.inlet[connection_index].V[i + 1] - flow.inlet[connection_index].V[i]).len();
  			}
  			flow.inlet[connection_index].l = l;
  		}
  		else {
  			flow.outlet[connection_index].V[0] = tmp_v1;
  			flow.outlet[connection_index].V[1] = tmp_v2;
  			// calculate new length
  			for (unsigned i = 0; i < flow.outlet[connection_index].V.size() - 1; i++) {
  				l += (flow.outlet[connection_index].V[i + 1] - flow.outlet[connection_index].V[i]).len();
  			}
  			flow.outlet[connection_index].l = l;
  		}
  
  		redisplay = true;
  		render_new_connection = false;
  		picked_connection = false;
  
  		flow.check_direct_connection();
  		flow.backup();
  	}

  }
  
  // define camera move based on mouse wheel move
  void glut_wheel(int wheel, int direction, int x, int y) {
  	

  	mods = glutGetModifiers();
  

  	mouse_x = x;
  	mouse_y = y;
  
  	GLdouble posX, posY, posZ;
  	window_to_world(posX, posY, posZ);			// get the world coordinates
  

  	if (!to_select_pressure && (simulation || build_inlet_outlet || manufacture)) {							// check current pixel position only in simualtion and build_inlet_outlet modes

  		bool flag = flow.epsilon_vertex((float)posX, (float)posY, (float)posZ, eps, scale, pressure_index);

  		if (flag && simulation && !glyph_mode) {

  			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);

  			undo = true;									// undo rendering

  		}

  		else if (!mods) {

  			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);
  		}
  	}

  
  	// rescale
  	if (mods == GLUT_ACTIVE_CTRL) {
  		if (direction > 0) {
  			if (scale >= 1)
  				scale += 1.0f;
  			else
  				scale += 0.1f;
  		}
  		else {
  			if (scale > 1)
  				scale -= 1.0f;
  			else if (scale <= 1 && scale > 0.1f)
  				scale -= 0.1f;
  			else
  				scale = 1.0f;
  		}

  		undo = true;
  		redisplay = true;

  	}

  	
  	glutPostRedisplay();
  }
  
  // define keyboard inputs
  void glut_keyboard(unsigned char key, int x, int y) {
  
  	// register different keyboard operation
  	switch (key) {

  		// saving network flow profile
  	case 's':
  		flow.save_network();

  		break;

  		// convert network to binary format (.nwt)
  	case 'c': {
  		std::vector<std::string> tmp = stim::parser::split(args.arg(0), '.');
  		std::stringstream ss;
  		ss << tmp[0] << ".nwt";
  		std::string filename = ss.str();
  		flow.saveNwt(filename);
  		break;
  	}

  
  		// subdivide current network for more detailed calculation
  	case 'd': {
  		// subdivide current network due to the limitation of current computation if needed
  		if (!subdivided && simulation && !glyph_mode) {
  			subdivided = true;
  
  			// check whether current network can be subdivided
  			if (flow.check_subdivision()) {
  				flow.subdivision();
  				get_background();
  			}
  
  			flow_initialize();					// initialize flow condition
  			// resimulation
  			flow_stable_state();				// main function of solving the linear system
  			flow.print_flow();
  			undo = true;
  		}
  		else if (subdivided && simulation && !glyph_mode) {
  			subdivided = false;
  			flow = backup;						// load back up
  
  			get_background();
  			flow_initialize();
  			// resimulation
  			flow_stable_state();				// main function of solving the linear system
  			flow.print_flow();
  			undo = true;
  		}
  		break;
  	}
  

  		// flow vector field visualization, Glyphs
  	case 'f':

  		if (glyph_mode && !manufacture && (simulation || build_inlet_outlet)) {
  			glyph_mode = false;
  			frame_mode = false;

  			redisplay = true;							// lines and arrows rendering use the same display list
  			int num = glutGet(GLUT_MENU_NUM_ITEMS);
  			if (num == 1)
  				glut_set_menu(num, 2);
  		}

  		else if (!glyph_mode && !manufacture && (simulation || build_inlet_outlet)) {
  			glyph_mode = true;

  			redisplay = true;
  			int num = glutGet(GLUT_MENU_NUM_ITEMS);
  			if (num == 2)
  				glut_set_menu(num, 1);
  		}

  		undo = true;
  		break;

  		// filaments around arrows
  	case 'g':
  		if (glyph_mode) {
  			if (frame_mode)
  				frame_mode = false;
  			else
  				frame_mode = true;
  		}

  		undo = true;

  		break;

  		// open/close index marks
  	case 'e':
  		if (mark_index)
  			mark_index = false;
  		else
  			mark_index = true;

  		undo = true;

  		break;
  

  		// output image stack
  	case 'm':
  		if (manufacture) {
  #ifdef __CUDACC__

  			flow.make_image_stack(S, dx, dy, dz, stackdir, image_stack, default_radius, scale);

  #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, in, out, 2, fragment_ratio);

  			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);

  	glutPassiveMotionFunc(glut_passive_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]);
  }
  

  // output an advertisement for the lab, authors and usage information
  void advertise() {
  	std::cout << std::endl << std::endl;
  	std::cout << " =======================================================================================" << std::endl;

  	std::cout << "|Thank you for using the AFAN tool!                                                     |" << std::endl;

  	std::cout << "|Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston                 |" << std::endl;
  	std::cout << "|Developers: Jiaming Guo, David Mayerich                                                |" << std::endl;
  	std::cout << "|Source: https://git.stim.ee.uh.edu/Jack/flow3.git									  |" << std::endl;
  	std::cout << " =======================================================================================" << std::endl << std::endl;
  

  	std::cout << "Usage(keyboard): e -> open/close indexing" << std::endl;
  	std::cout << "                 m -> build synthetic connections(connection mode)/output augmented network as image stack (manufacture mode)" << std::endl;

  	std::cout << "                 s -> save network flow profiles in profile folder as cvs files" << std::endl;
  	std::cout << "                 c -> convert .obj network to .nwt network and stores in main folder" << std::endl;
  	std::cout << "                 f -> open/close vector field visualization mode" << std::endl;
  	std::cout << "                 g -> render filament frames in vector fiedl visualization mode" << std::endl;

  	std::cout << args.str();
  }
  
  // main function: parse arguments and initialize GLUT

  int main(int argc, char* argv[]) {
  	
  	// add arguments
  	args.add("help", "prints the help");

  	args.add("units", "string indicating units of length for output measurements (ex. velocity)", "um", "text string");

  	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", "activate hilbert curves connections");

  	args.add("stack", "load the image stack");

  	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.add("scale", "scale down rendering fibers");
  	args.add("lcc", "extract the largest connected component");

  	
  	args.parse(argc, argv);								// parse the command line
  

  	if (args["help"].is_set()) {
  		advertise();
  		std::exit(1);
  	}
  

  	// load network

  	if (args.nargs() == 0) {
  		std::cout << "Network file required." << std::endl;
  		return 1;
  	}	
  	else {						// load network from user 
  		std::vector<std::string> tmp = stim::parser::split(args.arg(0), '.');

  		if ("obj" == tmp[1]) {

  			flow.load_obj(args.arg(0));

  			backup.load_obj(args.arg(0));
  		}	
  		else if ("nwt" == tmp[1]) {		// stim network binary format
  			flow.loadNwt(args.arg(0)); 
  			backup.loadNwt(args.arg(0));
  		}		
  		else if ("swc" == tmp[1]) {
  			flow.load_swc(args.arg(0)); 
  			backup.load_swc(args.arg(0));
  		}

  		else {
  			std::cout << "Invalid file type" << std::endl;
  			std::exit(1);
  		}
  	}

  
  	// extract the largest connected component

  	// get the units to work on
  	units = args["units"].as_string();
  	flow.set_units(units);
  

  	// blood pressure in capillaries range from 15 - 35 torr
  	// 1 torr = 133.3 Pa

  	max_pressure = (float)args["maxpress"].as_float();

  
  	// normal blood viscosity range from 4 - 15 mPa·s(cP)
  	// 1 Pa·s = 1 g / mm / s

  	u = (float)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 = (float)args["rou"].as_float();

  	// check whether to enable hilbert curves or not

  	hilbert_curve = args["hilbert"].is_set();

  	// load image stack if provided
  	if (args["stack"].is_set()) {
  		image_stack = true;
  		S.load_images(args["stack"].as_string());
  		// binary transformation
  #ifdef __CUDACC__

  		size_t N = S.samples();													// number of pixels loaded

  		unsigned char* d_S;														// image stack stored in device

  		unsigned char* h_S = (unsigned char*)malloc(N * sizeof(unsigned char));	// image stack stored in host
  		cudaMalloc((void**)&d_S, N * sizeof(unsigned char));
  		cudaMemcpy(d_S, S.data(), N * sizeof(unsigned char), cudaMemcpyHostToDevice);
  

  		size_t thread = 1024;
  		size_t block = N / thread + 1;

  		binary_transform <<<block, thread>>> (N, d_S, binary_threshold);		// binaryzation

  
  		cudaMemcpy(h_S, d_S, N * sizeof(unsigned char), cudaMemcpyDeviceToHost);
  
  		S.copy(h_S);

  #endif
  	}
  

  	// get the vexel and image stack size

  	dx = (float)args["stackres"].as_float(0);
  	dy = (float)args["stackres"].as_float(1);
  	dz = (float)args["stackres"].as_float(2);

  
  	// get the save directory of image stack
  	if (args["stackdir"].is_set())
  		stackdir = args["stackdir"].as_string();
  

  	// get the scale-down factor is provided
  	if (args["scale"].is_set())

  		scale = (float)args["scale"].as_float();

  	// glut main loop
  	bb = flow.boundingbox();
  	glut_initialize();
  	glutMainLoop();
  }