added a spherical harmonics visualization class, modified files to support it

David Mayerich
1 parent 6288483c
Showing 6 changed files with 364 additions and 69 deletions Show diff stats
stim/envi/binary.h
stim/envi/envi.h
stim/parser/arguments.h
stim/visualization/camera.h
stim/visualization/colormap.h
stim/visualization/sph_harmonics.h
@@ -225,6 +225,10 @@ public:
 		return true;
 	}
  
+	/// Reads a plane given a coordinate along the 0-axis (YZ plane)
+
+	/// @param p is a pointer to pre-allocated memory of size R[1] * R[2] * sizeof(T)
+	/// @param n is the 0-axis coordinate used to retrieve the plane
 	bool read_plane_0(T* p, unsigned int n){
  
 		if (n >= R[0]){										//make sure the number is within the possible range
@@ -247,6 +251,10 @@ public:
  
 	}
  
+	/// Reads a plane given a coordinate along the 1-axis (XZ plane)
+
+	/// @param p is a pointer to pre-allocated memory of size R[0] * R[2] * sizeof(T)
+	/// @param n is the 1-axis coordinate used to retrieve the plane
 	bool read_plane_1(T* p, unsigned int n){
  
 		unsigned int L = R[0] * sizeof(T);		//caculate the number of bytes in a sample line
@@ -268,10 +276,18 @@ public:
 		return true;
 	}
  
+	/// Reads a plane given a coordinate along the 2-axis (XY plane)
+
+	/// @param p is a pointer to pre-allocated memory of size R[0] * R[1] * sizeof(T)
+	/// @param n is the 2-axis coordinate used to retrieve the plane
 	bool read_plane_2(T* p, unsigned int n){
 		return read_page(p, n);
 	}
  
+	/// Reads a single pixel, treating the entire data set as a linear array
+
+	/// @param p is a pointer to pre-allocated memory of size sizeof(T)
+	/// @param i is the index to the pixel using linear indexing
 	bool read_pixel(T* p, unsigned int i){
 		if(i >= R[0] * R[1] * R[2]){
 			std::cout<<"ERROR read_pixel: n is out of range"<<std::endl;
@@ -283,6 +299,12 @@ public:
  
 	}
  
+	/// Reads a single pixel, given an x, y, z coordinate
+
+	/// @param p is a pointer to pre-allocated memory of size sizeof(T)
+	/// @param x is the x (0) axis coordinate
+	/// @param y is the y (1) axis coordinate
+	/// @param z is the z (2) axis coordinate
 	bool read_pixel(T* p, unsigned int x, unsigned int y, unsigned int z){
  
 		if(x < 0 || x >= R[0] || y < 0 || y >= R[1] || z < 0 || z > R[2]){
@@ -294,54 +316,6 @@ public:
 		return read_pixel(p, i);
 	}
  
-	//saves a hyperplane orthogonal to dimension d at intersection n
-	/*bool read_plane(T * dest, unsigned int d, unsigned int n){
-
-		//reset the file pointer back to the beginning of the file
-		file.seekg(0, std::ios::beg);
-
-		//compute the contiguous size C for each readable block
-		unsigned int C = 1;
-		for(unsigned int i = 0; i < d; i++)		//for each dimension less than d
-			C *= R[i];							//compute the product
-
-		//compute the non-contiguous size NC for each readable block
-		unsigned int NC = 1;
-		for(unsigned int i = d + 1; i < D; i++)
-			NC *= R[i];
-
-		//for all noncontiguous blocks, read each contiguous block that makes up the hyper-plane
-		for(unsigned int nc = 0; nc < NC; nc++){
-			file.seekg(n * C * sizeof(T), std::ios::cur);	//skip n contiguous blocks
-			file.read( (char*)&dest[nc * C], C * sizeof(T));	//read one contiguous block
-			file.seekg( (R[d] - n - 1) * C * sizeof(T), std::ios::cur);	//skip R[d] - n contiguous blocks
-		}
-
-		return true;
-
-	}*/
-
-	//save one pixel of the file into the memory, and return the pointer
-	/*bool read_spectrum(T * p, unsigned x, unsigned y){
-
-		unsigned int i;
-
-		if ( x >= R[0] || y >= R[1]){							//make sure the sample and line number is right
-			std::cout<<"ERROR: sample or line out of range"<<std::endl;
-			return false;
-		}
-
-		file.seekg((x + y * R[0]) * sizeof(T), std::ios::beg);           //point to the certain sample and line
-		for (i = 0; i < R[2]; i++)
-		{
-			file.read((char *)(p + i), sizeof(T));
-			file.seekg((R[1] * R[0] - 1) * sizeof(T), std::ios::cur);    //go to the next band
-		}
-
-		return true;	
-	}*/
-
-
 };
  
 }
@@ -710,6 +710,11 @@ public:
 		return false;
 	}
  
+	/// Retrieve a spectrum from the specified location
+
+	/// @param ptr is a pointer to pre-allocated memory of size B*sizeof(T)
+	/// @param x is the x-coordinate of the spectrum
+	/// @param y is the y-coordinate of the spectrum
 	bool spectrum(void* ptr, unsigned int x, unsigned int y){
  
 		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
@@ -256,14 +256,16 @@ namespace stim{
             args.add("foo", "foo takes a single integer value", "", "[intval]");
             args.add("bar", "bar takes two floating point values", "", "[value1], [value2]");
  
-        4) You generally want to immediately test for help and output available arguments:
+        4) Parse the command line:
+
+            args.parse(argc, argv);
+
+        5) You generally want to immediately test for help and output available arguments:
  
             if(args["help"].is_set())
                 std::cout<<args.str();
  
-        5) Parse the command line:
-
-            args.parse(argc, argv);
+        
  
         6)  Retrieve values:
  
@@ -436,7 +438,8 @@ namespace stim{
             }
  
             //set the last option
-            set(name, params);
+            if(name != "")
+            	set(name, params);
         }
  
         ///Determines of a parameter has been set and returns true if it has
@@ -167,14 +167,14 @@ public:
 	//output the camera settings
 	const void print(std::ostream& output)
 	{
-		output<<"Position: "<<p<<std::endl;
+		output<<"Position: "<<p.str()<<std::endl;
  
 	}
 	friend std::ostream& operator<<(std::ostream& out, const camera& c)
 	{
-		out<<"Position: "<<c.p<<std::endl;
-		out<<"Direction: "<<c.d<<std::endl;
-		out<<"Up: "<<c.up<<std::endl;
+		out<<"Position: "<<c.p.str()<<std::endl;
+		out<<"Direction: "<<c.d.str()<<std::endl;
+		out<<"Up: "<<c.up.str()<<std::endl;
 		out<<"Focal Distance: "<<c.focus<<std::endl;
 		return out;
 	}
@@ -221,7 +221,11 @@ static void cpuApplyBrewer(T* cpuSource, unsigned char* cpuDest, unsigned int N,
     for(int i=0; i<N; i++)
     {
         //compute the normalized value on [minVal maxVal]
-        float a = (cpuSource[i] - minVal) / (maxVal - minVal);
+		float a;
+		if(minVal != maxVal)
+			a = (cpuSource[i] - minVal) / (maxVal - minVal);
+		else
+			a = 0.5;
         if(a < 0) a = 0;
         if(a > 1) a = 1;
  
@@ -263,11 +267,15 @@ static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T 
         int i;
         float a;
         float range = valMax - valMin;
+
         for(i = 0; i<nVals; i++)
         {
             //normalize to the range [valMin valMax]
-            a = (cpuSource[i] - valMin) / range;
-
+			if(range != 0)
+				a = (cpuSource[i] - valMin) / range;
+			else
+				a = 0.5;
+	
             if(a < 0) a = 0;
             if(a > 1) a = 1;
  
@@ -279,22 +287,26 @@ static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T 
 }
  
 template<class T>
-static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, colormapType cm = cmGrayscale, bool positive = false)
+static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, colormapType cm = cmGrayscale)//, bool positive = false)
 {
     //computes the max and min range automatically
  
     //find the largest magnitude value
-    T maxVal = fabs(cpuSource[0]);
-    for(int i=0; i<nVals; i++)
+    T maxVal = cpuSource[0];
+    T minVal = cpuSource[0];
+    for(int i=1; i<nVals; i++)
 	{
-        if(fabs(cpuSource[i]) > maxVal)
-            maxVal = fabs(cpuSource[i]);
+        if(cpuSource[i] > maxVal)
+            maxVal = cpuSource[i];
+        if(cpuSource[i] < minVal)
+            minVal = cpuSource[i];
 	}
  
-    if(positive)
-        cpu2cpu(cpuSource, cpuDest, nVals, (T)0, maxVal, cm);
-    else
-        cpu2cpu(cpuSource, cpuDest, nVals, -maxVal, maxVal, cm);
+    //if(positive)
+    //    cpu2cpu(cpuSource, cpuDest, nVals, (T)0, maxVal, cm);
+    //else
+    //    cpu2cpu(cpuSource, cpuDest, nVals, -maxVal, maxVal, cm);
+    cpu2cpu(cpuSource, cpuDest, nVals, minVal, maxVal, cm);
  
 }
  
+#ifndef STIM_SPH_HARMONICS
+#define STIM_SPH_HARMONICS
+
+#include <GL/gl.h>
+
+#include <stim/gl/error.h>
+#include <stim/visualization/colormap.h>
+#include <vector>
+
+#define PI 3.14159
+#define WIRE_SCALE 1.001
+namespace stim{
+
+	class sph_harmonics{
+
+	private:
+
+		double* func;		//stores the raw function data (samples at each point)
+
+		GLuint color_tex;	//texture map that acts as a colormap for the spherical function
+		
+		unsigned int N;		//resolution of the spherical grid
+
+		std::vector<double> C;	//list of SH coefficients
+
+
+		//evaluates an associated Legendre polynomial (-l <= m <= l)
+		double P(int l,int m,double x)
+		{
+			// evaluate an Associated Legendre Polynomial P(l,m,x) at x
+			double pmm = 1.0;
+			if(m>0) {
+				double somx2 = sqrt((1.0-x)*(1.0+x));
+				double fact = 1.0;
+				for(int i=1; i<=m; i++) {
+					pmm *= (-fact) * somx2;
+					fact += 2.0;
+				}
+			}
+			if(l==m) return pmm;
+			double pmmp1 = x * (2.0*m+1.0) * pmm;
+			if(l==m+1) return pmmp1;
+			double pll = 0.0;
+			for(int ll=m+2; ll<=l; ++ll) {
+				pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m);
+				pmm = pmmp1;
+				pmmp1 = pll;
+			}
+			return pll;
+		}
+
+		//recursively calculate a factorial given a positive integer n
+		unsigned int factorial(unsigned int n) {
+			if (n == 0)
+			   return 1;
+			return n * factorial(n - 1);
+		}
+
+		//calculate the SH scaling constant
+		double K(int l, int m){
+
+			// renormalisation constant for SH function
+			double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*PI*factorial(l+m));
+			return sqrt(temp);
+		}
+
+		//calculate the value of the SH basis function (l, m) at (theta, phi)
+			//here, theta = [0, PI], phi = [0, 2*PI]
+		double SH(int l, int m, double theta, double phi){
+			// return a point sample of a Spherical Harmonic basis function
+			// l is the band, range [0..N]
+			// m in the range [-l..l]
+			// theta in the range [0..Pi]
+			// phi in the range [0..2*Pi]
+			const double sqrt2 = sqrt(2.0);
+			if(m==0) return K(l,0)*P(l,m,cos(theta));
+			else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
+			else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
+		}
+
+		void gen_function(){
+
+			//initialize the function to zero
+			memset(func, 0, sizeof(double) * N * N);
+
+			double theta, phi;
+			double result;
+			int l, m;
+
+			l = m = 0;
+			for(unsigned int c = 0; c < C.size(); c++){
+				
+
+				for(unsigned int xi = 0; xi < N; xi++)
+					for(unsigned int yi = 0; yi < N; yi++){
+
+						theta = (2 * PI) * ((double)xi / (N-1));
+						phi = PI * ((double)yi / (N-1));
+						result = C[c] * SH(l, m, phi, theta);		//phi and theta are reversed here (damn physicists)
+						func[yi * N + xi] += result;
+					}
+
+				m++;			//increment m
+
+				//if we're in a new tier, increment l and set m = -l
+				if(m > l){		
+					l++;
+					m = -l;
+				}
+			}
+		}
+
+		void gl_prep_draw(){
+
+			//enable depth testing
+				//this has to be used instead of culling because the sphere can have negative values
+			glEnable(GL_DEPTH_TEST);
+			glDepthMask(GL_TRUE);
+			glEnable(GL_TEXTURE_2D);	//enable 2D texture mapping
+		}
+
+		//draw a texture mapped sphere representing the function surface
+		void gl_draw_sphere() {
+
+			//PI is used to convert from spherical to cartesian coordinates
+			//const double PI = 3.14159;
+
+			//bind the 2D texture representing the color map
+			glBindTexture(GL_TEXTURE_2D, color_tex);
+
+			//Draw the Sphere
+			int i, j;
+
+			for(i = 1; i <= N-1; i++) {
+				double phi0 = PI * ((double) (i - 1) / (N-1));
+				double phi1 = PI * ((double) i / (N-1));
+
+				glBegin(GL_QUAD_STRIP);
+				for(j = 0; j <= N; j++) {
+
+					//calculate the indices into the function array
+					int phi0_i = i-1;
+					int phi1_i = i;
+					int theta_i = j;
+					if(theta_i == N)
+						theta_i = 0;
+
+					double v0 = func[phi0_i * N + theta_i];
+					double v1 = func[phi1_i * N + theta_i];
+
+					v0 = fabs(v0);
+					v1 = fabs(v1);
+
+
+					double theta = 2 * PI * (double) (j - 1) / N;
+					double x0 = v0 * cos(theta) * sin(phi0);
+					double y0 = v0 * sin(theta) * sin(phi0);
+					double z0 = v0 * cos(phi0);
+
+					double x1 = v1 * cos(theta) * sin(phi1);
+					double y1 = v1 * sin(theta) * sin(phi1);
+					double z1 = v1 * cos(phi1);
+
+					glTexCoord2f(theta / (2 * PI), phi0 / PI);
+					glVertex3f(x0, y0, z0);
+
+					glTexCoord2f(theta / (2 * PI), phi1 / PI);
+					glVertex3f(x1, y1, z1);
+				}
+				glEnd();
+			}
+		}
+
+		//draw a wire frame sphere representing the function surface
+		void gl_draw_wireframe() {
+
+					//PI is used to convert from spherical to cartesian coordinates
+					//const double PI = 3.14159;
+
+					//bind the 2D texture representing the color map
+					glDisable(GL_TEXTURE_2D);
+					glColor3f(0.0f, 0.0f, 0.0f);
+
+					//Draw the Sphere
+					int i, j;
+
+					for(i = 1; i <= N-1; i++) {
+						double phi0 = PI * ((double) (i - 1) / (N-1));
+						double phi1 = PI * ((double) i / (N-1));
+
+						glBegin(GL_LINE_STRIP);
+						for(j = 0; j <= N; j++) {
+
+							//calculate the indices into the function array
+							int phi0_i = i-1;
+							int phi1_i = i;
+							int theta_i = j;
+							if(theta_i == N)
+								theta_i = 0;
+
+							double v0 = func[phi0_i * N + theta_i];
+							double v1 = func[phi1_i * N + theta_i];
+
+							v0 = fabs(v0);
+							v1 = fabs(v1);
+
+
+							double theta = 2 * PI * (double) (j - 1) / N;
+							double x0 = WIRE_SCALE * v0 * cos(theta) * sin(phi0);
+							double y0 = WIRE_SCALE * v0 * sin(theta) * sin(phi0);
+							double z0 = WIRE_SCALE * v0 * cos(phi0);
+
+							double x1 = WIRE_SCALE * v1 * cos(theta) * sin(phi1);
+							double y1 = WIRE_SCALE * v1 * sin(theta) * sin(phi1);
+							double z1 = WIRE_SCALE * v1 * cos(phi1);
+
+							glTexCoord2f(theta / (2 * PI), phi0 / PI);
+							glVertex3f(x0, y0, z0);
+
+							glTexCoord2f(theta / (2 * PI), phi1 / PI);
+							glVertex3f(x1, y1, z1);
+						}
+						glEnd();
+					}
+				}
+
+		void init(unsigned int n){
+
+			//set the sphere resolution
+			N = n;
+
+			//allocate space for the color map
+			unsigned int bytes = N * N * sizeof(unsigned char) * 3;
+			unsigned char* color_image;
+			color_image = (unsigned char*) malloc(bytes);
+
+			//allocate space for the function
+			func = (double*) malloc(N * N * sizeof(double));
+
+			//generate a function (temporary)
+			gen_function();
+
+			//generate a colormap from the function
+			stim::cpu2cpu<double>(func, color_image, N*N, stim::cmBrewer);
+
+			//prep everything for drawing
+			gl_prep_draw();			
+
+			//generate an OpenGL texture map in the current context
+			glGenTextures(1, &color_tex);
+			//bind the texture
+			glBindTexture(GL_TEXTURE_2D, color_tex);
+
+			//copy the color data from the buffer to the GPU
+			glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, color_image);
+
+			//initialize all of the texture parameters
+			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
+			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_REPLACE);
+			
+			//free the buffer
+			free(color_image);
+		}
+
+
+	public:
+
+		void glRender(){
+			//set all OpenGL parameters required for drawing
+			gl_prep_draw();
+
+			//draw the sphere
+			gl_draw_sphere();
+			//gl_draw_wireframe();
+
+		}
+
+		void glInit(unsigned int n){
+			init(n);
+		}
+
+		void push(double c){
+			C.push_back(c);
+		}
+
+
+
+
+
+	};		//end class sph_harmonics
+
+
+
+
+}
+
+
+#endif