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::gl_spharmonics<double> S;

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)

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	
	S.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
	S.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 = 2*stim::PI;	///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");

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

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

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



}