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;



}