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.01
  #define phi_scale		0.01

  #define zoom_scale		0.1

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

  stim::spharmonics<double> S;
  stim::gl_spharmonics<double> R;				//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 zaxis = 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.0;					//alpha value for interpolation
  float dalpha = 0.1;					//change in alpha value with a key press
  

  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;
  }
  
  //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]);
  
  	//draw the sphere

  	if (interp)
  		R = (S * (1.0f - alpha) + Si * alpha);
  	else
  		R = S;
  	R.glRender();
  

  	//glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT);
  
  	//draw the z-axis if requested
  	if(zaxis){
  		glDisable(GL_TEXTURE_2D);
  		glColor3f(0.0f, 1.0f, 0.0f);
  		glBegin(GL_LINES);
  			glVertex3f(0.0, 0.0, 0.0);
  			glVertex3f(0.0, 0.0, 100.0);
  		glEnd();
  	}
  

  	//flush commands on the GPU
  	glutSwapBuffers();
  }
  
  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 the mouse is pressed
  	if(state == GLUT_DOWN){
  		//set the current mouse position
  		mx = x;		my = y;
  	}
  }
  
  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 = 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] = 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(time(NULL));			///set random seed
  	srand(100);			///set random seed
  
  	for(int i = 0; i < num_samples; i++)
  	{
  		z = uniformRandom()*range + Z[1];

  		theta = uniformRandom() * 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;
  }
  

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

  	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.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("zaxis", "render the z-axis as a green line");

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

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

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

  	//set the z-axis flag
  	if(args["zaxis"].is_set())
  		zaxis = true;
  

  	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.nargs() > 0){
  		//push all of the arguments to the spherical harmonics class as coefficients
  		for(unsigned int a = 0; a < args.nargs(); a++)
  			S.push(atof(args.arg(a).c_str()));
  	}
  
  	//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(time(NULL));
  
  		unsigned int N = args["rand"].as_int(0);		//get the number of random coefficients
  		double Cmin = args["rand"].as_float(1);			//get the minimum and maximum coefficient values
  		double Cmax = args["rand"].as_float(2);
  
  		//generate the coefficients
  		for(unsigned int c = 0; c < N; c++){
  
  			double norm = (double) rand() / RAND_MAX;		//calculate a random number in the range [0, 1]
  			double scaled = norm * (Cmax - Cmin) + Cmin;	//scale the random number to [Cmin, Cmax]
  			S.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
  		S.resize(Cmax + 1);
  		
  		//insert each coefficient
  		for(unsigned int i = 0; i < nC; i++){
  			S.setc(C[i], V[i]);
  		}
  
  	}
  	else if(args["obj"].is_set()){
  
  		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
  	else 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++)
  			S.push(0);
  
  		//add the n'th coefficient
  		S.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();
  		}
  	}
  
  
  	
  	

  
  	//if the user asks for help, give it and exit
  	if(args["help"].is_set()){
  		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);
  	}
  }
  
  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");
  
  	//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;
  
  
  
  }