stim/visualization/gl_spharmonics.h

#ifndef STIM_GL_SPHARMONICS_H
#define STIM_GL_SPHARMONICS_H

#include <stim/math/spharmonics.h>
#include <stim/gl/error.h>
#include <stim/math/vec3.h>
#include <stim/math/constants.h>
#include <stim/visualization/colormap.h>

namespace stim {

	template<typename T>
	class gl_spharmonics {

		GLuint dlist;
		GLuint tex;

		bool displacement;
		bool colormap;
		bool magnitude;

		void init_tex() {
			T* sfunc = (T*)malloc(N * N * sizeof(T));								//create a 2D array to store the spherical function
			Sc.get_func(sfunc, N, N);														//generate the spherical function based on the Sc coefficients
			unsigned char* tex_buffer = (unsigned char*)malloc(3 * N * N);			//create a buffer to store the texture map
			stim::cpu2cpu<T>(sfunc, tex_buffer, N * N, stim::cmBrewer);				//create a Brewer colormap based on the spherical function
			stim::buffer2image(tex_buffer, "sfunc.ppm", N, N);

			if (tex) glDeleteTextures(1, &tex);															//if a texture already exists, delete it
			glGenTextures(1, &tex);																		//create a new texture and store the ID																
			glBindTexture(GL_TEXTURE_2D, tex);															//bind the texture
			glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, tex_buffer);		//copy the color data from the buffer to the GPU
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);								//initialize all of the texture parameters
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
			glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
		}

	public:
		stim::spharmonics<T> Sc;					//spherical harmonic representing the color component
		stim::spharmonics<T> Sd;					//spherical harmonic representing the displacement component
		size_t N;

		gl_spharmonics(size_t slices) {
			N = slices;
			dlist = 0;								//initialize the display list index to zero (no list)
			tex = 0;								//initialize the texture index to zero (no texture)
			displacement = true;
			colormap = true;
			magnitude = true;
		}

		gl_spharmonics(stim::spharmonics<T> disp, stim::spharmonics<T> color, size_t slices) :
			 gl_spharmonics<T>(slices)
		{
			Sc = color;
			Sd = disp;
		}

		~gl_spharmonics() {
			if (dlist) glDeleteLists(dlist, 1);		//delete the display list when the object is destroyed
		}

                ///saves the coeffients of a spherical harmonics in a text documents.
                void save_coeffs(std::string filename)
                {
                        std::ofstream myfile;
                        myfile.open(filename.c_str());
                        for(int i = 0; i < Sd.C.size()-1; i++)
                        {
                                myfile << Sd.C[i] << ",";
                        }
                        myfile << Sd.C[Sd.C.size()-1] << std::endl;
                        for(int i = 0; i < Sc.C.size()-1; i++)
                        {
                                myfile << Sc.C[i] << ",";
                        }
                        myfile << Sc.C[Sc.C.size()-1] << std::endl;
                        myfile.close();
                }

                void save(std::string filename)
                {

                        if (!tex) {
                                init_tex();
                        }

                        stim::obj<T> object;
                        size_t theta_i, phi_i;
                        T d_theta = (T)stim::TAU / (T)N;
                        T d_phi = (T)stim::PI / (T)(N-1);
                        object.matKd("sfunc.jpg");
                        for (phi_i = 1; phi_i < N; phi_i++) {
                                T phi = phi_i * d_phi;
                                object.Begin(OBJ_TRIANGLE_STRIP);
                                for (theta_i = 0; theta_i <= N; theta_i++) {
                                        T theta = (N - theta_i) * d_theta;
                                        float theta_t = 1 - (float)theta_i / (float)N;

                                        T r;
                                        if (!displacement) r = 1;                                               //if no displacement, set the r value to 1 (renders a unit sphere)
                                        else r = Sd(theta, phi);                                                //otherwise calculate the displacement value

                                        glColor3f(1.0f, 1.0f, 1.0f);
                                        if (!colormap) {                                                                                //if no colormap is being rendered
                                                if (r < 0) glColor3f(1.0, 0.0, 0.0);                            //if r is negative, render it red
                                                else glColor3f(0.0, 1.0, 0.0);                                          //otherwise render in green
                                        }
                                        if (magnitude) {                                                                                //if the magnitude is being displaced, calculate the magnitude of r
                                                if (r < 0) r = -r;
                                        }
                                        stim::vec3<T> s(r, theta, phi);
                                        stim::vec3<T> c = s.sph2cart();
                                        stim::vec3<T> n;                                                                                                                //allocate a value to store the normal
                                        if (!displacement) n = c;                                                                                               //if there is no displacement, the normal is spherical
                                        else n = Sd.dphi(theta, phi).cross(Sd.dtheta(theta, phi));                              //otherwise calculate the normal as the cross product of derivatives

                                        object.TexCoord(theta_t, (float)phi_i / (float)N);
                                        //std::cout << theta_t <<"        "<<(float)phi_i / (float)N << "----------------";
                                        object.Normal(n[0], n[1], n[2]);
                                        object.Vertex(c[0], c[1], c[2]);

                                        T r1;
                                        if (!displacement) r1 = 1;
                                        else r1 = Sd(theta, phi - d_phi);
                                        if (!colormap) {                                                                                //if no colormap is being rendered
                                                if (r1 < 0)     glColor3f(1.0, 0.0, 0.0);                               //if r1 is negative, render it red
                                                else glColor3f(0.0, 1.0, 0.0);                                          //otherwise render in green
                                        }
                                        if (magnitude) {                                                                                //if the magnitude is being rendered, calculate the magnitude of r
                                                if (r1 < 0) r1 = -r1;
                                        }
                                        stim::vec3<T> s1(r1, theta, phi - d_phi);
                                        stim::vec3<T> c1 = s1.sph2cart();
                                        stim::vec3<T> n1;
                                        if (!displacement) n1 = c1;
                                        else n1 = Sd.dphi(theta, phi - d_phi).cross(Sd.dtheta(theta, phi - d_phi));

                                        //std::cout << theta_t << "        " << (float)(phi_i - 1) / (float)N << std::endl;
                                        object.TexCoord(theta_t, 1.0/(2*(N)) + (float)(phi_i-1) / (float)N);
                                        object.Normal(n1[0], n1[1], n1[2]);
                                        object.Vertex(c1[0], c1[1], c1[2]);
                                }
                                object.End();
                        }
                        object.matKd();
                        object.save(filename);
                }


		/// This function renders the spherical harmonic to the current OpenGL context
		void render() {
			//glShadeModel(GL_FLAT);
			glPushAttrib(GL_ENABLE_BIT);
			glDisable(GL_CULL_FACE);
			glEnable(GL_DEPTH_TEST);
			glDepthMask(GL_TRUE);
			
			if (!tex) {
				init_tex();
			}

			if (colormap) {
				glEnable(GL_TEXTURE_2D);
				glBindTexture(GL_TEXTURE_2D, tex);
			}

			if (!dlist) {
				dlist = glGenLists(1);
				glNewList(dlist, GL_COMPILE);
				glPushAttrib(GL_ENABLE_BIT);
				glEnable(GL_NORMALIZE);

				//Draw the Sphere
				size_t theta_i, phi_i;
				T d_theta = (T)stim::TAU / (T)N;
				T d_phi = (T)stim::PI / (T)(N-1);

				for (phi_i = 1; phi_i < N; phi_i++) {
					T phi = phi_i * d_phi;
					glBegin(GL_QUAD_STRIP);
					for (theta_i = 0; theta_i <= N; theta_i++) {
						T theta = (N - theta_i) * d_theta;
						float theta_t = 1 - (float)theta_i / (float)N;

						T r;
						if (!displacement) r = 1;						//if no displacement, set the r value to 1 (renders a unit sphere)
						else r = Sd(theta, phi);						//otherwise calculate the displacement value

						glColor3f(1.0f, 1.0f, 1.0f);
						if (!colormap) {										//if no colormap is being rendered
							if (r < 0) glColor3f(1.0, 0.0, 0.0);				//if r is negative, render it red
							else glColor3f(0.0, 1.0, 0.0);						//otherwise render in green
						}
						if (magnitude) {										//if the magnitude is being displaced, calculate the magnitude of r
							if (r < 0) r = -r;
						}
						stim::vec3<T> s(r, theta, phi);
						stim::vec3<T> c = s.sph2cart();
						stim::vec3<T> n;														//allocate a value to store the normal
						if (!displacement) n = c;												//if there is no displacement, the normal is spherical
						else n = Sd.dphi(theta, phi).cross(Sd.dtheta(theta, phi));				//otherwise calculate the normal as the cross product of derivatives

						glTexCoord2f(theta_t, (float)phi_i / (float)N);
						//std::cout << theta_t <<"        "<<(float)phi_i / (float)N << "----------------";
						glNormal3f(n[0], n[1], n[2]);
						glVertex3f(c[0], c[1], c[2]);

						T r1;
						if (!displacement) r1 = 1;
						else r1 = Sd(theta, phi - d_phi);
						if (!colormap) {										//if no colormap is being rendered
							if (r1 < 0)	glColor3f(1.0, 0.0, 0.0);				//if r1 is negative, render it red
							else glColor3f(0.0, 1.0, 0.0);						//otherwise render in green
						}
						if (magnitude) {										//if the magnitude is being rendered, calculate the magnitude of r
							if (r1 < 0) r1 = -r1;
						}
						stim::vec3<T> s1(r1, theta, phi - d_phi);
						stim::vec3<T> c1 = s1.sph2cart();					
						stim::vec3<T> n1;
						if (!displacement) n1 = c1;
						else n1 = Sd.dphi(theta, phi - d_phi).cross(Sd.dtheta(theta, phi - d_phi));

						//std::cout << theta_t << "        " << (float)(phi_i - 1) / (float)N << std::endl;
						glTexCoord2f(theta_t, 1.0/(2*(N)) + (float)(phi_i-1) / (float)N);
						glNormal3f(n1[0], n1[1], n1[2]);
						glVertex3f(c1[0], c1[1], c1[2]);
					}
					glEnd();
				}
				glPopAttrib();
				glEndList();
			}
			glCallList(dlist);											//call the display list to render

			glPopAttrib();
			glDisable(GL_DEPTH_TEST);
			glDepthMask(GL_FALSE);
		}

		/// Push a coefficient to the spherical harmonic - by default, push applies the component to both the displacement and color SH
		void push(T coeff) {
			Sd.push(coeff);
			Sc.push(coeff);
		}

		/// Resize the spherical harmonic coefficient array
		void resize(size_t s) {
			Sd.resize(s);
			Sc.resize(s);
		}

		/// Set a spherical harmonic coefficient to the given value
		void setc(unsigned int c, T value) {
			Sd.setc(c, value);
			Sc.setc(c, value);
		}

		void project(T* data, size_t x, size_t y, size_t nc) {
			Sd.project(data, x, y, nc);
			Sc = Sd;
		}

		/// Project a set of samples onto the basis
		void project(std::vector<vec3<float> >& vlist, size_t nc) {
			Sd.project(vlist, nc);
			Sc = Sd;
		}

		/// Calculate a density function from a list of points in spherical coordinates
		void pdf(std::vector<stim::vec3<T> >& vlist, size_t nc) {
			Sd.pdf(vlist, nc);
			Sc = Sd;
		}

		void slices(size_t s) {
			N = s;
		}

		size_t slices() {
			return N;
		}

		void rendermode(bool displace, bool color, bool mag = true) {
			displacement = displace;
			colormap = color;
			magnitude = mag;
		}

	};
}


#endif