diff --git a/stim/envi/binary.h b/stim/envi/binary.h index 5925cea..61bd3b1 100644 --- a/stim/envi/binary.h +++ b/stim/envi/binary.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"<= 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"< 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 1) a = 1; @@ -279,22 +287,26 @@ static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T } template -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 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); } diff --git a/stim/visualization/sph_harmonics.h b/stim/visualization/sph_harmonics.h new file mode 100644 index 0000000..bc0a540 --- /dev/null +++ b/stim/visualization/sph_harmonics.h @@ -0,0 +1,301 @@ +#ifndef STIM_SPH_HARMONICS +#define STIM_SPH_HARMONICS + +#include + +#include +#include +#include + +#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 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(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 -- libgit2 0.21.4