visualization / shview

  #include <iostream>

  #include <sstream>

  
  #include <GL/glut.h>
  
  #include <stim/visualization/camera.h>
  #include <stim/parser/arguments.h>

  #include <stim/visualization/obj.h>
  #include <stim/visualization/gl_spharmonics.h>

  #include <stim/math/constants.h>

  #define theta_scale		0.01f
  #define phi_scale		0.01f
  #define zoom_scale		0.1f

  
  //create a global camera that will specify the viewport
  stim::camera cam;
  int mx, my;			//mouse coordinates in the window space
  

  stim::gl_spharmonics<float> SH(300);				//spherical harmonics to render

  
  float d = 1.5;		//initial distance between the camera and the sphere
  
  bool rotate_zoom = true;	//sets the current camera mode (rotation = true, zoom = false)
  
  stim::arglist args;			//class for processing command line arguments
  

  bool axis = false;			//render the z-axis (set via a command line flag)

  /// INTERPOLATION

  //stim::gl_spharmonics<double> Si;

  bool interp = false;				//if we are interpolating

  float alpha = 0.0f;					//alpha value for interpolation
  float dalpha = 0.1f;					//change in alpha value with a key press
  
  /// Visualization flags
  bool mag = true;
  bool cmap = true;
  bool displace = true;
  bool light = true;

  bool init(){
  
  	//set the clear color to white
  	glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
  
  	//initialize the camera
  	cam.setPosition(d, d, d);
  	cam.LookAt(0, 0, 0, 0, 1, 1);
  	cam.setFOV(40);
  
  	//initialize the texture map stuff	

  	//R.glInit(256);

  
  	return true;
  }
  

  void render_axes() {
  	glDisable(GL_TEXTURE_2D);							//turn off texture mapping		
  	glBegin(GL_LINES);
  		glColor3f(1.0f, 0.0f, 0.0f);					//set the color to RED and render X
  		glVertex3f(0.0, 0.0, 0.0);
  		glVertex3f(100.0, 0.0, 0.0);
  
  		glColor3f(0.0f, 1.0f, 0.0f);					//set the color to RED and render X
  		glVertex3f(0.0, 0.0, 0.0);
  		glVertex3f(0.0, 100.0, 0.0);
  
  		glColor3f(0.0f, 0.0f, 1.0f);					//set the color to RED and render X
  		glVertex3f(0.0, 0.0, 0.0);
  		glVertex3f(0.0, 0.0, 100.0);
  	glEnd();
  }
  
  void render_lights() {
  
  	stim::vec3<float> view = cam.getDirection();							//get the camera view direction
  	stim::vec3<float> up = cam.getUp();										//get the camera up direction
  	stim::vec3<float> left = view.cross(up).norm();							//create a vector pointing to the left
  	
  	GLfloat light0_diffuse[] = { 0.5f, 0.5f, 0.5f, 1 };						//create a directional light shining from the viewer's left
  	GLfloat light0_position[] = { left[0], left[1], left[2], 0.0 };
  	glLightfv(GL_LIGHT0, GL_DIFFUSE, light0_diffuse);
  	glLightfv(GL_LIGHT0, GL_POSITION, light0_position);
  
  	GLfloat light1_diffuse[] = { 0.8f, 0.8f, 0.8f, 1 };						//create a directional light shining from the viewer's position
  	GLfloat light1_position[] = { -view[0], -view[1], -view[2], 0.0 };
  	glLightfv(GL_LIGHT1, GL_DIFFUSE, light0_diffuse);
  	glLightfv(GL_LIGHT1, GL_POSITION, light1_position);
  	
  
  }
  

  //code that is run every time the user changes something
  void display(){
  	//clear the screen
  	glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT);
  
  	//set the projection matrix
  	glMatrixMode(GL_PROJECTION);				//put the projection matrix on the stack
  	glLoadIdentity();							//set it to the identity matrix
  	gluPerspective(cam.getFOV(), 1, 0.001, 1000000);	//set up a perspective projection
  
  
  	//set the model view matrix
  	glMatrixMode(GL_MODELVIEW);					//load the model view matrix to the stack
  	glLoadIdentity();							//set it to the identity matrix
  
  	//get the camera parameters

  	stim::vec3<float> p = cam.getPosition();
  	stim::vec3<float> u = cam.getUp();
  	stim::vec3<float> d = cam.getDirection();

  
  	//specify the camera parameters to OpenGL
  	gluLookAt(p[0], p[1], p[2], d[0], d[1], d[2], u[0], u[1], u[2]);
  

  	if (light) {
  		glEnable(GL_LIGHTING);
  		glEnable(GL_LIGHT0);
  		glEnable(GL_LIGHT1);
  	}
  	
  	glEnable(GL_COLOR_MATERIAL);
  	glEnable(GL_CULL_FACE);
  	glEnable(GL_DEPTH_TEST);
  
  	render_lights();									//render the lights
  	glColor3f(1.0, 1.0, 1.0);
  	SH.render();
  
  	if (light) {
  		glDisable(GL_LIGHTING);								//disable lighting to render the axes
  		glDisable(GL_LIGHT0);
  		glDisable(GL_LIGHT1);

  	}
  

  	if(axis) render_axes();								//render the axes if the user requests them
  	glutSwapBuffers();									//swap in the back buffer (double-buffering is used to prevent tearing)

  }
  
  void mouse_press(int button, int state, int x, int y){
  
  	//set the camera motion mode based on the mouse button pressed
  	if(button == GLUT_LEFT_BUTTON)
  		rotate_zoom = true;

  	else if(button == GLUT_RIGHT_BUTTON)

  		rotate_zoom = false;

  	if(state == GLUT_DOWN){								//if the mouse is pressed
  		mx = x;		my = y;								//set the current mouse position

  	}
  }
  
  void mouse_drag(int x, int y){
  
  	//if the camera is in rotation mode, rotate
  	if(rotate_zoom == true){
  		float theta = theta_scale * (mx - x);
  		float phi = -phi_scale * (my - y);
  
  		//if the mouse is dragged
  		cam.OrbitFocus(theta, phi);
  	}
  	//otherwize zoom
  	else{

  		cam.Push(zoom_scale*(my - y));

  	}
  
  	//update the mouse position
  	mx = x;		my = y;
  
  	glutPostRedisplay();
  }
  

  float uniformRandom()
  {
  	return (  (float)(rand()))/(  (float)(RAND_MAX)); 
  }
  
  std::vector<stim::vec3 <float> >

  sample_sphere(int num_samples, float radius = 1.0) {

  	float solidAngle = (float)stim::TAU;	///Solid angle to sample over
  	float PHI[2], Z[2], range;		///Range of angles in cylinderical coordinates
  	PHI[0] = solidAngle/2;			///project the solid angle into spherical coords
  	PHI[1] = (float)asin(0);		
  	Z[0] = cos(PHI[0]);				///project the z into spherical coordinates
  	Z[1] = cos(PHI[1]);		
  	range = Z[0] - Z[1];			///the range of all possible z values.

  	float z, theta, phi;			/// temporary individual

  
  	std::vector<stim::vec3<float> > samples;
  

  	srand(100);						///set random seed

  
  	for(int i = 0; i < num_samples; i++)
  	{
  		z = uniformRandom()*range + Z[1];

  		theta = uniformRandom() * (float)stim::TAU;

  		phi = acos(z);
  		stim::vec3<float> sph(1, theta, phi);
  		stim::vec3<float> cart = sph.sph2cart();
  		sph[0] *= radius;
  		samples.push_back(cart);
  	}
  	samples.push_back(stim::vec3<float>(0.,0.,1.));
  	samples.push_back(stim::vec3<float>(0.,1.0,0.));
  	samples.push_back(stim::vec3<float>(0.,-1.,0.));
  
  
  	std::stringstream name;
        for(int i = 0; i < num_samples; i++)
             name << samples[i].str() << std::endl;
             name << samples[num_samples].str() << std::endl;
             name << samples[num_samples+1].str() << std::endl;
             name << samples[num_samples+2].str() << std::endl;
  	
      
        std::ofstream outFile;
        outFile.open("New_Pos_Vectors.txt");
        outFile << name.str().c_str();
  
  	return samples;
  }
  

  std::vector<stim::vec3<float>> obj2vec3(stim::obj<float>& o) {
  	stim::vec3<float> c = o.centroid();									//calculate the centroid
  	size_t nv = o.numV();												//get the number of vertices in the obj file
  	std::vector<stim::vec3<float>> vlist(nv);							//create an array to store the vertices
  	stim::vec3<float> p, v, s;
  	for (size_t vi = 0; vi < nv; vi++) {								//for each vertex in the obj file
  		v = o.getV(vi);
  		p = (v - c);													//center and normalize
  		s = p.cart2sph();												//convert to spherical coordinates
  		vlist[vi] = s;													//store the spherical coordinates in the point list
  	}
  	return vlist;
  }
  
  void advertise() {
  	std::cout << "usage: shview c0 c1 c2 c3 ... --option [A B C]" << std::endl;
  	std::cout << "examples:" << std::endl;
  	std::cout << "   generate a spherical function with 4 coefficients (l=0 to 2)" << std::endl;
  	std::cout << "          shview 1.3 0.2 2.3 1.34" << std::endl;
  	std::cout << "   display a spherical function representing the spherical harmonic l = 3, m = -2" << std::endl;
  	std::cout << "          shview --basis 3 -2" << std::endl;
  	std::cout << args.str();
  	exit(0);
  }
  

  void process_arguments(int argc, char* argv[]){
  
  	args.add("help", "prints this help");

  	args.section("Coefficients");
  	args.add("coef", "specify spherical harmonic coefficients", "", "c0 c1 c2 c3 ...");

  	args.add("rand", "generates a random set of SH coefficients", "", "[N min max]");
  	args.add("sparse", "generates a function based on a set of sparse basis functions", "", "[l0 m0 c0 l1 m1 c1 l2 m2 c2 ...]");

  	args.section("Projection");
  	args.add("image", "approximates a spherical function represented by an image", "", "filename1.ppm filename2.ppm");
  	args.add("n", "number of spherical harmonics coefficients to use for a projection", "10", "positive integer");	
  	args.add("basis", "displays the specified SH basis function", "", "n, or [l m]");	

  	args.add("obj", "approximates a geometric object given as a Wavefront OBJ file", "", "filename");
  	args.add("out", "filename for outputting spherical harmonics coefficients", "", "filename");

  	args.add("pdf", "outputs the PDF if an OBJ files is given");

  	args.section("Visualization");
  	args.add("axis", "render the z-axis as a green line");
  	args.add("nodisp", "render the spherical function without displacement");
  	args.add("nocmap", "render the spherical function without color mapping");
  	args.add("nomag", "render the spherical function without calculating the absolute value (negative values are displaced negatively");
  	args.add("nolight", "render the spherical function without lighting");

  	args.add("interp", "interpolates between two specified sets of coefficients", "", "[c0 c1 c2 c3 ...]");

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

  	//if the user asks for help, give it and exit
  	if (args["help"].is_set()) {
  		advertise();

  	}
  

  	//set the z-axis flag
  	if(args["axis"].is_set())
  		axis = true;
  	if (args["nodisp"]) displace = false;
  	if (args["nocmap"]) cmap = false;
  	if (args["nomag"]) mag = false;
  	if (args["nolight"]) light = false;
  	SH.rendermode(displace, cmap, mag);
  
  	//if (args["interp"]) {											//if an interpolation harmonic is provided
  	//	for (unsigned int a = 0; a < args["interp"].nargs(); a++)
  	//		Si.push(args["interp"].as_float(a));
  	//}
  

  	//if arguments are specified, push them as coefficients

  	if(args["coef"]){

  		//push all of the arguments to the spherical harmonics class as coefficients

  		for(unsigned int a = 0; a < args["coef"].nargs(); a++)
  			SH.push((float)args["coef"].as_float(a));

  	}

  	//if the user wants to use a random set of SH coefficients
  	else if(args["rand"].is_set()){
  
  		//return an error if the user specifies both fixed and random coefficients
  		if(args.nargs() != 0){
  			std::cout<<"Error: both fixed and random coefficients are specified"<<std::endl;
  			exit(1);
  		}
  
  		//seed the random number generator

  		srand((unsigned int)time(NULL));

  		unsigned int N = args["rand"].as_int(0);				//get the number of random coefficients
  		float Cmin = (float)args["rand"].as_float(1);			//get the minimum and maximum coefficient values
  		float Cmax = (float)args["rand"].as_float(2);

  
  		//generate the coefficients
  		for(unsigned int c = 0; c < N; c++){
  

  			float norm = (float) rand() / RAND_MAX;		//calculate a random number in the range [0, 1]
  			float scaled = norm * (Cmax - Cmin) + Cmin;	//scale the random number to [Cmin, Cmax]
  			SH.push(scaled);									//push the value as a coefficient

  		}
  	}

  	else if (args["sparse"].is_set()) {

  
  		//calculate the number of sparse coefficients
  		unsigned int nC = args["sparse"].nargs() / 3;
  
  		std::vector<unsigned int> C;	//vector of 1D coefficients
  		unsigned int Cmax = 0;			//maximum coefficient provided
  
  		std::vector<double> V;			//vector of 1D coefficient values
  
  		unsigned int c;
  		int l, m;
  		double v;
  		//for each provided coefficient

  		for (unsigned int i = 0; i < nC; i++) {

  
  			//load data for a single coefficient from the command line

  			l = args["sparse"].as_int(i * 3 + 0);
  			m = args["sparse"].as_int(i * 3 + 1);
  			v = args["sparse"].as_float(i * 3 + 2);

  
  			//calculate the 1D coefficient
  			c = pow(l + 1, 2) - (l - m) - 1;
  
  			//update the maximum coefficient index

  			if (c > Cmax) Cmax = c;

  			//insert the coefficient and value into vectors
  			C.push_back(c);

  			V.push_back(v);

  		}
  
  		//set the size of the SH coefficient array

  		SH.resize(Cmax + 1);
  

  		//insert each coefficient

  		for (unsigned int i = 0; i < nC; i++) {
  			SH.setc(C[i], V[i]);
  		}
  	}
  	else if (args["image"]) {
  		if (args["image"].nargs() == 0) {
  			std::cout << "ERROR: an image file must be specified as an argument to --image" << std::endl;
  			exit(1);
  		}
  		stim::image<float> I(args["image"].as_string(0));							//load the image
  		if (I.channels() > 1) I = I.channel(0);										//if a color image is provided, convert to monochrome
  		I = I.stretch(0.0f, 1.0f);													//stretch the function so that it lies in [0 1]
  		SH.project(I.data(), I.width(), I.height(), args["n"].as_int(0));			//project the image function onto a spherical harmonics basis
  		if (args["image"].nargs() > 1) {											//if the user provides two images
  			I.load(args["image"].as_string(1));
  			if (I.channels() > 1) I = I.channel(0);									//if a color image is provided, convert to monochrome
  			I = I.stretch(0.0f, 1.0f);												//stretch the function so that it lies in [0 1]
  			SH.Sc.project(I.data(), I.width(), I.height(), args["n"].as_int(0));

  		}

  	}

  	else if(args["obj"].is_set()){

  		if (args["obj"].nargs() == 0) {
  			std::cout << "ERROR: an OBJ file must be specified as an argument to --obj" << std::endl;
  			exit(1);
  		}
  		stim::obj<float> OBJ(args["obj"].as_string(0));								//open the OBJ file
  		std::vector<stim::vec3<float>> vlist = obj2vec3(OBJ);				//get the centered points for the OBJ
  		SH.pdf(vlist, args["n"].as_int(0));		

  		/*std::string filename = args["obj"].as_string(0);

  		unsigned int l = args["obj"].as_int(1);

  		int p = args["obj"].as_int(2);
  		std::cout << p << std::endl;
  		std::vector<stim::vec3<float> > sphere = sample_sphere(p);
  		p = p+3;

  
  		//create an obj object
  		stim::obj<double> object(filename);
  
  		//get the centroid of the object
  		stim::vec<double> c = object.centroid();

  		c[0] = 0; c[1] = 0; c[2] = 0;

  
  		//get the number of vertices in the model
  		unsigned int nV = object.numV();
  
  		//for each vertex in the model, create an MC sample
  		std::vector< stim::vec<double> > spherical;
  		stim::vec<float> sample;
  		stim::vec<float> centered;
  		for(unsigned int i = 0; i < nV; i++){
  
  			sample = object.getV(i);			//get a vertex in cartesian coordinates
  			centered = sample - c;
  			spherical.push_back(centered.cart2sph());
  		}
  
  		//generate the spherical PDF
  		stim::spharmonics<double> P;
  		P.pdf(spherical, l, l);

  		std::vector<float> weights;		///array of weights
  		if(args["pdf"].is_set())
  		{
  //			S.pdf(spherical, l, l);
  			for(int i = 0; i < p; i++)	///for each point on the sphere.
  			{
  				float val = 0;		///value starts with 0
  				for(int j = 0; j < nV; j++)		///for each point on surface
  				{
  					stim::vec3<float> star(object.getV(j)[0] - c[0],
  						object.getV(j)[1] - c[1], 
  						object.getV(j)[2] - c[2]);		///center each point on the model 
  //					val += abs(star.dot(sphere[i])); 		///sum the dot product of the centered point and the sphere.
  					if(star.dot(sphere[i]) > 0)
  						val += pow(star.dot(sphere[i]),8); 		///sum the dot product of the centered point and the sphere.
  				}
  				weights.push_back(val);
  			}
  			
  			S.mcBegin(l,l);
  			for(int i = 0; i < p; i++)
  			{
  				if(sphere[i] == stim::vec3<float>(0., 0., 1.))
  				{
  					std::cout << i  << sphere[i] << " " << weights[i] << std::endl;
  				}
  				if(sphere[i] == stim::vec3<float>(0., 1., 0.))
  				{
  					std::cout << i << sphere[i] << " " << weights[i] << std::endl;
  				}
  				if(sphere[i] == stim::vec3<float>(0., -1., 0.))
  				{
  					std::cout << i <<  sphere[i] << " " << weights[i] << std::endl;
  				}
  				stim::vec3<float> sph = sphere[i].cart2sph();
  				S.mcSample(sph[1], sph[2], weights[i]);
  			}
  			S.mcEnd();

  			*/

  		}

  		/*else{ 

  			//begin Monte-Carlo sampling, using the model vertices as samples
  			S.mcBegin(l, l);
  			double theta, phi, fx, px;
  			for(unsigned int i = 0; i < nV; i++){
  				theta = spherical[i][1];
  				phi = spherical[i][2];
  				fx = spherical[i][0];
  				px = P(theta, phi);
  				S.mcSample(theta, phi, fx / px);
  			}
  			S.mcEnd();

  		}*/
  	//}

  
  	//if the user specifies an SH basis function

  	/*if(args["basis"].is_set()){

  
  		unsigned int n;
  
  		//if the user specifies one index for the basis function
  		if(args["basis"].nargs() == 1)
  			n = args["basis"].as_int(0);
  		else if(args["basis"].nargs() == 2){
  			int l = args["basis"].as_int(0);		//2D indexing (l, m)
  			int m = args["basis"].as_int(1);
  
  			n = pow(l+1, 2) - (l - m) - 1;			//calculate the 1D index
  		}
  
  		//add zeros for the first (n-1) coefficients
  		for(unsigned int c = 0; c < n; c++)

  			SH.push(0);

  
  		//add the n'th coefficient

  		SH.push(1);
  	}*/

  
  	//output the spherical harmonics coefficients if requested

  	/*if(args["out"].is_set()){

  		
  		if(args["out"].nargs() == 0)
  			std::cout<<S.str()<<std::endl;
  		else{
  
  			//open the output file
  			std::ofstream outfile;
  			outfile.open(args["out"].as_string(0).c_str());
  
  			outfile<<S.str();
  
  			outfile.close();
  		}

  	}*/

  }
  
  int main(int argc, char *argv[]){
  

  #ifdef _WIN32
  	args.set_ansi(false);
  #endif
  

  	//initialize GLUT
  	glutInit(&argc, argv);
  
  	//process arguments
  	process_arguments(argc, argv);
  
  	//set the size of the GLUT window
  	glutInitWindowSize(500, 500);
  
  	glutInitDisplayMode(GLUT_DEPTH | GLUT_RGBA | GLUT_DOUBLE);
  
  	//create the GLUT window (and an OpenGL context)
  	glutCreateWindow("Spherical Harmonic Viewport");
  

  	//SH.rendermode(false, true, true);
  	//SH.colormap();										//use color-mapping
  

  	//set the display function (which will be called repeatedly by glutMainLoop)
  	glutDisplayFunc(display);
  
  	//set the mouse press function (called when a mouse button is pressed)
  	glutMouseFunc(mouse_press);
  	//set the mouse motion function (which will be called any time the mouse is dragged)
  	glutMotionFunc(mouse_drag);
  
  	//run the initialization function
  	if(!init())
  		return 1;	//return an error if it fails
  
  
  	//enter the main loop
  	glutMainLoop();
  
  	//return 0 if everything is awesome
  	return 0;
  
  
  
  }