From 01707489309c7fd07f3372e43063de936141e4cb Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Sun, 12 Feb 2017 23:12:37 -0600 Subject: [PATCH] implemented a new spherical harmonics renderer --- stim/image/image.h | 13 ++++++------- stim/math/spharmonics.h | 486 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- stim/visualization/gl_spharmonics-dep.h | 259 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ stim/visualization/gl_spharmonics.h | 374 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 4 files changed, 711 insertions(+), 421 deletions(-) create mode 100644 stim/visualization/gl_spharmonics-dep.h diff --git a/stim/image/image.h b/stim/image/image.h index f686e48..36482a8 100644 --- a/stim/image/image.h +++ b/stim/image/image.h @@ -137,10 +137,6 @@ public: std::cout << "Error opening input file in image::load_netpbm()" << std::endl; exit(1); } - if (sizeof(T) != 1) { - std::cout << "Error in image::load_netpbm() - data type must be 8-bit integer." << std::endl; - exit(1); - } size_t nc; //allocate space for the number of channels char format[2]; //allocate space to hold the image format tag @@ -187,9 +183,12 @@ public: } size_t h = atoi(sh.c_str()); //convert the string into an integer - allocate(w, h, nc); //allocate space for the image - infile.read((char*)img, size()); //copy the binary data from the file to the image - infile.close(); + allocate(w, h, nc); //allocate space for the image + unsigned char* buffer = (unsigned char*)malloc(w * h * nc); //create a buffer to store the read data + infile.read((char*)buffer, size()); //copy the binary data from the file to the image + infile.close(); //close the file + for (size_t n = 0; n < size(); n++) img[n] = (T)buffer[n]; //copy the buffer data into the image + free(buffer); //free the buffer array } diff --git a/stim/math/spharmonics.h b/stim/math/spharmonics.h index e6c5a11..a2fc711 100644 --- a/stim/math/spharmonics.h +++ b/stim/math/spharmonics.h @@ -2,280 +2,348 @@ #define STIM_SPH_HARMONICS #include -#include #include #include #include #include #define WIRE_SCALE 1.001 -namespace stim{ +namespace stim { -template -class spharmonics{ + template + class spharmonics { -public: - std::vector C; //list of SH coefficients + public: + std::vector C; //list of SH coefficients -protected: + protected: + unsigned int mcN; //number of Monte-Carlo samples + unsigned int coeff_1d(unsigned int l, int m) { //convert (l,m) to i (1D coefficient index) + return pow(l + 1, 2) - (l - m) - 1; + } + void coeff_2d(size_t c, unsigned int& l, int& m) { //convert a 1D coefficient index into (l, m) + l = (unsigned int)ceil(sqrt((double)c + 1)) - 1; //the major index is equal to sqrt(c) - 1 + m = (int)(c - (size_t)(l * l)) - (int)l; //the minor index is calculated by finding the difference + } + public: + spharmonics() { + mcN = 0; + } + spharmonics(size_t c) : spharmonics() { + resize(c); + } - unsigned int mcN; //number of Monte-Carlo samples + void push(T c) { + C.push_back(c); + } - //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){ - //std::complex result = boost::math::spherical_harmonic(l, m, phi, theta); - //return result.imag() + result.real(); - - //this calculation is based on calculating the real spherical harmonics: - // https://en.wikipedia.org/wiki/Spherical_harmonics#Addition_theorem - if (m < 0) { - return sqrt(2.0) * pow(-1, m) * boost::math::spherical_harmonic(l, abs(m), phi, theta).imag(); + void resize(unsigned int n) { + C.resize(n); } - else if (m == 0) { - return boost::math::spherical_harmonic(l, m, phi, theta).real(); + + void setc(unsigned int l, int m, T value) { + unsigned int c = coeff_1d(l, m); + C[c] = value; } - else { - return sqrt(2.0) * pow(-1, m) * boost::math::spherical_harmonic(l, m, phi, theta).real(); + + T getc(unsigned int l, int m) { + unsigned int c = coeff_1d(l, m); + return C[c]; } - } - unsigned int coeff_1d(unsigned int l, int m){ - return pow(l + 1, 2) - (l - m) - 1; - } + void setc(unsigned int c, T value) { + C[c] = value; + } - + unsigned int getSize() const { + return C.size(); + } + std::vector getC() const { + return C; + } + //calculate the value of the SH basis function (l, m) at (theta, phi) + //here, theta = [0, PI], phi = [0, 2*PI] + T SH(unsigned int l, int m, T theta, T phi) { + //std::complex result = boost::math::spherical_harmonic(l, m, phi, theta); + //return result.imag() + result.real(); + + //this calculation is based on calculating the real spherical harmonics: + // https://en.wikipedia.org/wiki/Spherical_harmonics#Addition_theorem + if (m < 0) { + return sqrt(2.0) * pow(-1, m) * boost::math::spherical_harmonic(l, abs(m), phi, theta).imag(); + } + else if (m == 0) { + return boost::math::spherical_harmonic(l, m, phi, theta).real(); + } + else { + return sqrt(2.0) * pow(-1, m) * boost::math::spherical_harmonic(l, m, phi, theta).real(); + } + } -public: - spharmonics() { - mcN = 0; - } - spharmonics(size_t c) : spharmonics() { - resize(c); - } + /// Calculate the spherical harmonic result given a 1D coefficient index + T SH(size_t c, T theta, T phi) { + unsigned int l; + int m; + coeff_2d(c, l, m); + return SH(l, m, theta, phi); + } - void push(double c){ - C.push_back(c); - } - void resize(unsigned int n){ - C.resize(n); - } - void setc(unsigned int l, int m, T value){ - unsigned int c = coeff_1d(l, m); - C[c] = value; - } + /// Initialize Monte-Carlo sampling of a function using N spherical harmonics coefficients - void setc(unsigned int c, T value){ - C[c] = value; - } + /// @param N is the number of spherical harmonics coefficients used to represent the user function + void mcBegin(unsigned int coefficients) { + C.resize(coefficients, 0); + mcN = 0; + } - unsigned int getSize() const{ - return C.size(); - } + void mcBegin(unsigned int l, int m) { + unsigned int c = pow(l + 1, 2) - (l - m); + mcBegin(c); + } - std::vector getC() const{ - return C; - } + void mcSample(T theta, T phi, T val) { - + int l, m; + T sh; - /// Initialize Monte-Carlo sampling of a function using N spherical harmonics coefficients + l = m = 0; + for (unsigned int i = 0; i < C.size(); i++) { - /// @param N is the number of spherical harmonics coefficients used to represent the user function - void mcBegin(unsigned int coefficients){ - C.resize(coefficients, 0); - mcN = 0; - } + sh = SH(l, m, theta, phi); + C[i] += sh * val; - void mcBegin(unsigned int l, int m){ - unsigned int c = pow(l + 1, 2) - (l - m); - mcBegin(c); - } + m++; //increment m - void mcSample(double theta, double phi, double val){ + //if we're in a new tier, increment l and set m = -l + if (m > l) { + l++; + m = -l; + } + } //end for all coefficients - int l, m; - double sh; + //increment the number of samples + mcN++; - l = m = 0; - for(unsigned int i = 0; i < C.size(); i++){ + } //end mcSample() - sh = SH(l, m, theta, phi); - C[i] += sh * val; + void mcEnd() { - m++; //increment m + //divide all coefficients by the number of samples + for (unsigned int i = 0; i < C.size(); i++) + C[i] /= mcN; + } - //if we're in a new tier, increment l and set m = -l - if(m > l){ - l++; - m = -l; + /// Generates a PDF describing the probability distribution of points on a spherical surface + /// @param sph_pts is a list of points in spherical coordinates (theta, phi) where theta = [0, 2pi] and phi = [0, pi] + /// @param l is the maximum degree of the spherical harmonic function + /// @param m is the maximum order + void pdf(std::vector > sph_pts, unsigned int l, int m) { + mcBegin(l, m); //begin spherical harmonic sampling + unsigned int nP = sph_pts.size(); + for (unsigned int p = 0; p < nP; p++) { + mcSample(sph_pts[p][1], sph_pts[p][2], 1.0); } - } //end for all coefficients + mcEnd(); + } + + void pdf(std::vector > sph_pts, size_t c) { + unsigned int l; + int m; + coeff_2d(c, l, m); + pdf(sph_pts, l, m); + } + + /// Project a set of samples onto a spherical harmonic basis + void project(std::vector > sph_pts, unsigned int l, int m) { + mcBegin(l, m); //begin spherical harmonic sampling + unsigned int nP = sph_pts.size(); + for (unsigned int p = 0; p < nP; p++) { + mcSample(sph_pts[p][1], sph_pts[p][2], sph_pts[p][0]); + } + mcEnd(); + } + void project(std::vector > sph_pts, size_t c) { + unsigned int l; + int m; + coeff_2d(c, l, m); + project(sph_pts, l, m); + } - //increment the number of samples - mcN++; + /// Generates a PDF describing the density distribution of points on a sphere + /// @param sph_pts is a list of points in cartesian coordinates + /// @param l is the maximum degree of the spherical harmonic function + /// @param m is the maximum order + /// @param c is the centroid of the points in sph_pts. DEFAULT 0,0,0 + /// @param n is the number of points of the surface of the sphere used to create the PDF. DEFAULT 1000 + /// @param norm, a boolean that sets where the output vectors will be normalized between 0 and 1. + /// @param + /*void pdf(std::vector > sph_pts, unsigned int l, int m, stim::vec3 c = stim::vec3(0, 0, 0), unsigned int n = 1000, bool norm = true, std::vector w = std::vector()) + { + std::vector weights; ///the weight at each point on the surface of the sphere. + // weights.resize(n); + unsigned int nP = sph_pts.size(); + std::vector > sphere = stim::Random::sample_sphere(n, 1.0, stim::TAU); + if (w.size() < nP) + w = std::vector(nP, 1.0); + + for (int i = 0; i < n; i++) + { + T val = 0; + for (int j = 0; j < nP; j++) + { + stim::vec3 temp = sph_pts[j] - c; + if (temp.dot(sphere[i]) > 0) + val += pow(temp.dot(sphere[i]), 4)*w[j]; + } + weights.push_back(val); + } - } //end mcSample() + mcBegin(l, m); //begin spherical harmonic sampling - void mcEnd(){ + if (norm) + { + T min = *std::min_element(weights.begin(), weights.end()); + T max = *std::max_element(weights.begin(), weights.end()); + for (unsigned int i = 0; i < n; i++) + { + stim::vec3 sph = sphere[i].cart2sph(); + mcSample(sph[1], sph[2], (weights[i] - min) / (max - min)); + } - //divide all coefficients by the number of samples - for(unsigned int i = 0; i < C.size(); i++) - C[i] /= mcN; - } + } + else { + for (unsigned int i = 0; i < n; i++) + { + stim::vec3 sph = sphere[i].cart2sph(); + mcSample(sph[1], sph[2], weights[i]); + } + } + mcEnd(); + }*/ - /// Generates a PDF describing the probability distribution of points on a spherical surface - /// @param sph_pts is a list of points in spherical coordinates (theta, phi) where theta = [0, 2pi] and phi = [0, pi] - /// @param l is the maximum degree of the spherical harmonic function - /// @param m is the maximum order + std::string str() { - void pdf(std::vector > sph_pts, unsigned int l, int m){ - - mcBegin( l, m ); //begin spherical harmonic sampling + std::stringstream ss; - unsigned int nP = sph_pts.size(); + int l, m; + l = m = 0; + for (unsigned int i = 0; i < C.size(); i++) { - for(unsigned int p = 0; p < nP; p++){ - mcSample(sph_pts[p][1], sph_pts[p][2], 1.0); - } + ss << C[i] << '\t'; - mcEnd(); - } - - /// Generates a PDF describing the probability distribution of points on a spherical surface - /// @param sph_pts is a list of points in cartesian coordinates - /// @param l is the maximum degree of the spherical harmonic function - /// @param m is the maximum order - /// @param c is the centroid of the points in sph_pts. DEFAULT 0,0,0 - /// @param n is the number of points of the surface of the sphere used to create the PDF. DEFAULT 1000 - void pdf(std::vector > sph_pts, unsigned int l, int m, stim::vec3 c = stim::vec3(0,0,0), unsigned int n = 1000) - { - std::vector weights; ///the weight at each point on the surface of the sphere. -// weights.resize(n); - unsigned int nP = sph_pts.size(); - std::vector > sphere = stim::Random::sample_sphere(n, 1.0, stim::TAU); - for(int i = 0; i < n; i++) - { - double val = 0; - for(int j = 0; j < nP; j++) - { - stim::vec3 temp = sph_pts[j] - c; - if(temp.dot(sphere[i]) > 0) - val += pow(temp.dot(sphere[i]),4); - } - weights.push_back(val); - } - - - mcBegin(l, m); //begin spherical harmonic sampling - - for(unsigned int i = 0; i < n; i++) - { - stim::vec3 sph = sphere[i].cart2sph(); - mcSample(sph[1], sph[2], weights[i]); - } + m++; //increment m - mcEnd(); - } + //if we're in a new tier, increment l and set m = -l + if (m > l) { + l++; + m = -l; - std::string str(){ + ss << std::endl; - std::stringstream ss; + } + } - int l, m; - l = m = 0; - for(unsigned int i = 0; i < C.size(); i++){ - - ss< l){ - l++; - m = -l; + } - ss< l) { + l++; + m = -l; + } } + return fx; } - return ss.str(); + /// Returns the derivative of the spherical function with respect to theta + /// return value is in cartesian coordinates + vec3 dtheta(T theta, T phi, T d = 0.01) { + T r = p(theta, phi); //calculate the value of the spherical function at three points + T rt = p(theta + d, phi); + //double rp = p(theta, phi + d); + vec3 s(r, theta, phi); //get the spherical coordinate position for all three points + vec3 st(rt, theta + d, phi); + //vec3 sp(rp, theta, phi + d); - } + vec3 c = s.sph2cart(); + vec3 ct = st.sph2cart(); + //vec3 cp = sp.sph2cart(); - /// Returns the value of the function at the coordinate (theta, phi) + vec3 dt = (ct - c)/d; //calculate the derivative + return dt; + } - /// @param theta = [0, 2pi] - /// @param phi = [0, pi] - double operator()(double theta, double phi){ + /// Returns the derivative of the spherical function with respect to phi + /// return value is in cartesian coordinates + vec3 dphi(T theta, T phi, T d = 0.01) { + T r = p(theta, phi); //calculate the value of the spherical function at three points + //double rt = p(theta + d, phi); + T rp = p(theta, phi + d); - double fx = 0; + vec3 s(r, theta, phi); //get the spherical coordinate position for all three points + //vec3 st(rt, theta + d, phi); + vec3 sp(rp, theta, phi + d); - int l = 0; - int m = 0; - for(unsigned int i = 0; i < C.size(); i++){ - fx += C[i] * SH(l, m, theta, phi); - m++; - if(m > l){ - l++; - m = -l; - } + vec3 c = s.sph2cart(); + //vec3 ct = st.sph2cart(); + vec3 cp = sp.sph2cart(); + vec3 dp = (cp - c) / d; //calculate the derivative + return dp; } - - return fx; - } - - /// Fill an NxN grid with the spherical function for theta = [0 2pi] and phi = [0 pi] - void grid(T* data, size_t N){ - double dt = stim::TAU / (double)N; //calculate the step size in each direction - double dp = stim::PI / (double)(N - 1); - for(size_t ti = 0; ti < N; ti++){ - for(size_t pi = 0; pi < N){ - data[pi * N + ti] = (*this)((double)ti * dt, (double)pi * dp); - } + + /// Returns the value of the function at the coordinate (theta, phi) + /// @param theta = [0, 2pi] + /// @param phi = [0, pi] + T operator()(T theta, T phi) { + return p(theta, phi); } - } -/* - //overload arithmetic operations - spharmonics operator*(T rhs) const { - spharmonics result(C.size()); //create a new spherical harmonics object - for (size_t c = 0; c < C.size(); c++) //for each coefficient - result.C[c] = C[c] * rhs; //calculate the factor and store the result in the new spharmonics object - return result; - } - - spharmonics operator+(spharmonics rhs) { - size_t low = std::min(C.size(), rhs.C.size()); //store the number of coefficients in the lowest object - size_t high = std::max(C.size(), rhs.C.size()); //store the number of coefficients in the result - bool rhs_lowest = false; //true if rhs has the lowest number of coefficients - if (rhs.C.size() < C.size()) rhs_lowest = true; //if rhs has a lower number of coefficients, set the flag - - spharmonics result(high); //create a new object - size_t c; - for (c = 0; c < low; c++) //perform the first batch of additions - result.C[c] = C[c] + rhs.C[c]; //perform the addition - - for (c = low; c < high; c++) { - if (rhs_lowest) - result.C[c] = C[c]; - else - result.C[c] = rhs.C[c]; + /// Fill an NxN grid with the spherical function for theta = [0 2pi] and phi = [0 pi] + void get_func(T* data, size_t X, size_t Y) { + T dt = stim::TAU / (T)X; //calculate the step size in each direction + T dp = stim::PI / (T)(Y - 1); + for (size_t ti = 0; ti < X; ti++) { + for (size_t pi = 0; pi < Y; pi++) { + data[pi * X + ti] = (*this)((T)ti * dt, (T)pi * dp); + } + } } - return result; - } - spharmonics operator-(spharmonics rhs) { - return (*this) + (rhs * (T)(-1)); - } -*/ -}; //end class sph_harmonics + /// Project a spherical function onto the basis using C coefficients + /// @param data is a pointer to the function values in (theta, phi) coordinates + /// @param N is the number of samples along each axis, where theta = [0 2pi), phi = [0 pi] + void project(T* data, size_t x, size_t y, size_t nc) { + stim::cpu2image(data, "test.ppm", x, y, stim::cmBrewer); + C.resize(nc, 0); //resize the coefficient array to store the necessary coefficients + T dtheta = stim::TAU / (T)(x - 1); //calculate the grid spacing along theta + T dphi = stim::PI / (T)y; //calculate the grid spacing along phi + T theta, phi; + for (size_t c = 0; c < nc; c++) { //for each coefficient + for (size_t theta_i = 0; theta_i < x; theta_i++) { //for each coordinate in the provided array + theta = theta_i * dtheta; //calculate theta + for (size_t phi_i = 0; phi_i < y; phi_i++) { + phi = phi_i * dphi; //calculate phi + C[c] += data[phi_i * x + theta_i] * SH(c, theta, phi) * dtheta * dphi * sin(phi); + } + } + } + } + }; //end class sph_harmonics @@ -283,4 +351,4 @@ public: } -#endif +#endif \ No newline at end of file diff --git a/stim/visualization/gl_spharmonics-dep.h b/stim/visualization/gl_spharmonics-dep.h new file mode 100644 index 0000000..37f2414 --- /dev/null +++ b/stim/visualization/gl_spharmonics-dep.h @@ -0,0 +1,259 @@ +#ifndef STIM_GL_SPHARMONICS_H +#define STIM_GL_SPHARMONICS_H + +#include + +#include +#include +#include + +namespace stim{ + + +template +class gl_spharmonics : public spharmonics{ + +protected: + + using stim::spharmonics::SH; + stim::spharmonics surface; //harmonic that stores the surface information + stim::spharmonics color; //harmonic that stores the color information + T* func_surface; //stores the raw function data (samples at each point) + T* func_color; //stores the raw color 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 + + ///generates the + void gen_surface(){ + //allocate memory and initialize the function to zero + func_surface = (T*) malloc(N * N * sizeof(T)); + memset(func_surface, 0, sizeof(T) * N * N); + + double theta, phi; + double result; + int l, m; + + + l = m = 0; + //for each coefficient + for(unsigned int c = 0; c < surface.C.size(); c++){ + + //iterate through the entire 2D grid representing the function + for(unsigned int xi = 0; xi < N; xi++){ + for(unsigned int yi = 0; yi < N; yi++){ + + //get the spherical coordinates for each grid point + theta = (2 * PI) * ((double)xi / (N-1)); + phi = PI * ((double)yi / (N-1)); + + //sum the contribution of the current spherical harmonic based on the coefficient + result = surface.C[c] * SH(l, m, theta, phi); + + //store the result in a 2D array (which will later be used as a texture map) + func_surface[yi * N + xi] += result; + } + } + + //keep track of m and l here + m++; //increment m + //if we're in a new tier, increment l and set m = -l + if(m > l){ + l++; + m = -l; + } + } + } + + ///generates the color function + void gen_color(){ + + //initialize the function to zero + func_color = (T*) malloc(N * N * sizeof(T)); + memset(func_color, 0, sizeof(T) * N * N); + + double theta, phi; + double result; + int l, m; + + + l = m = 0; + //for each coefficient + for(unsigned int c = 0; c < color.C.size(); c++){ + + //iterate through the entire 2D grid representing the function + for(unsigned int xi = 0; xi < N; xi++){ + for(unsigned int yi = 0; yi < N; yi++){ + + //get the spherical coordinates for each grid point + theta = (2 * PI) * ((double)xi / (N-1)); + phi = PI * ((double)yi / (N-1)); + + //sum the contribution of the current spherical harmonic based on the coefficient + result = color.C[c] * SH(l, m, theta, phi); + + //store the result in a 2D array (which will later be used as a texture map) + func_color[yi * N + xi] += result; + } + } + + //keep track of m and l here + 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() { + + //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_surface[phi0_i * N + theta_i]; + double v1 = func_surface[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(); + } + + glBindTexture(GL_TEXTURE_2D, 0); + } + + void gen_texture() + { + //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); + + //generate a colormap from the function + stim::cpu2cpu(func_color, 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); + } + + void gl_init(unsigned int n){ + + //set the sphere resolution + N = n; + + //set up the surface data. + gen_surface(); + + //set up the color data and generate the corresponding texture. + gen_color(); + gen_texture(); + + } + +public: + gl_spharmonics(unsigned int n = 256) { + gl_init(n); + } + + gl_spharmonics( stim::spharmonics in_function, unsigned int n = 256) + { + surface = in_function; + color = in_function; + gl_init(n); + } + + gl_spharmonics( stim::spharmonics in_function, stim::spharmonics in_color, unsigned int n = 256) + { + surface = in_function; + color = in_color; + gl_init(n); + } + + gl_spharmonics& operator=(const spharmonics rhs) { + gl_spharmonics result(rhs.C.size()); + result.C = spharmonics::rhs.C; + return result; + } + + void glRender(){ + //set all OpenGL parameters required for drawing + gl_prep_draw(); + + //draw the sphere + gl_draw_sphere(); + } + + void glInit(unsigned int n = 256){ + gl_init(n); + } +}; //end gl_spharmonics + + +}; //end namespace stim + + + + +#endif diff --git a/stim/visualization/gl_spharmonics.h b/stim/visualization/gl_spharmonics.h index 921cc50..9882325 100644 --- a/stim/visualization/gl_spharmonics.h +++ b/stim/visualization/gl_spharmonics.h @@ -1,235 +1,199 @@ #ifndef STIM_GL_SPHARMONICS_H #define STIM_GL_SPHARMONICS_H -#include - +#include #include +#include +#include #include -#include - -namespace stim{ - - -template -class gl_spharmonics : public spharmonics{ - -protected: - using spharmonics::C; - using spharmonics::SH; - - T* 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 - - void gen_function(){ - - //initialize the function to zero - memset(func, 0, sizeof(double) * N * N); - - double theta, phi; - double result; - int l, m; +namespace stim { + + template + 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(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); + } - l = m = 0; - //for each coefficient - for(unsigned int c = 0; c < C.size(); c++){ + public: + stim::spharmonics Sc; //spherical harmonic representing the color component + stim::spharmonics 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; + } - //iterate through the entire 2D grid representing the function - for(unsigned int xi = 0; xi < N; xi++){ - for(unsigned int yi = 0; yi < N; yi++){ + ~gl_spharmonics() { + if (dlist) glDeleteLists(dlist, 1); //delete the display list when the object is destroyed + } - //get the spherical coordinates for each grid point - theta = (2 * PI) * ((double)xi / (N-1)); - phi = PI * ((double)yi / (N-1)); + /// This function renders the spherical harmonic to the current OpenGL context + void render() { + //glShadeModel(GL_FLAT); + glPushAttrib(GL_ENABLE_BIT); + glDisable(GL_CULL_FACE); + + if (!tex) { + init_tex(); + } - //sum the contribution of the current spherical harmonic based on the coefficient - result = C[c] * SH(l, m, theta, phi); + if (colormap) { + glEnable(GL_TEXTURE_2D); + glBindTexture(GL_TEXTURE_2D, tex); + } - //store the result in a 2D array (which will later be used as a texture map) - func[yi * N + xi] += result; + 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 s(r, theta, phi); + stim::vec3 c = s.sph2cart(); + stim::vec3 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 s1(r1, theta, phi - d_phi); + stim::vec3 c1 = s1.sph2cart(); + stim::vec3 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 - //keep track of m and l here - m++; //increment m - //if we're in a new tier, increment l and set m = -l - if(m > l){ - l++; - m = -l; - } + glPopAttrib(); } - } - - 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() { - - //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(); + /// 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); } - } - - void gl_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 functional representation that will be used for the texture map and vertices - 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); + /// Resize the spherical harmonic coefficient array + void resize(size_t s) { + Sd.resize(s); + Sc.resize(s); + } - //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: - gl_spharmonics() { - - } - - gl_spharmonics(const spharmonics copy) : gl_spharmonics() { - C = copy.C; - } - - gl_spharmonics& operator=(const spharmonics rhs) { - gl_spharmonics result(rhs.C.size()); - result.C = spharmonics::rhs.C; - return result; - } - - void glRender(){ - //set all OpenGL parameters required for drawing - gl_prep_draw(); - - //draw the sphere - gl_draw_sphere(); - } - - void glInit(unsigned int n = 256){ - gl_init(n); - gen_function(); - } - - //overload arithmetic operations - gl_spharmonics operator*(T rhs) const { - gl_spharmonics result(C.size()); //create a new sph erical harmonics object - for (size_t c = 0; c < C.size(); c++) //for each coefficient - result.C[c] = C[c] * rhs; //calculat e the factor and store the result in the new spharmonics object - return result; - } - - gl_spharmonics operator+(gl_spharmonics rhs) { - size_t low = std::min(C.size(), rhs.C.size()); //store th e number of coefficients in the lowest object - size_t high = std::max(C.size(), rhs.C.size()); //store th e number of coefficients in the result - bool rhs_lowest = false; //true if rhs has the lowest number of coefficients - if (rhs.C.size() < C.size()) rhs_lowest = true; //if rhs has a low er number of coefficients, set the flag - - gl_spharmonics result(high); //create a new object - size_t c; - for (c = 0; c < low; c++) //perform the first batch of additions - result.C[c] = C[c] + rhs.C[c]; // perform the addition - - for (c = low; c < high; c++) { - if (rhs_lowest) - result.C[c] = C[c]; - else - result.C[c] = rhs.C[c]; + /// Set a spherical harmonic coefficient to the given value + void setc(unsigned int c, T value) { + Sd.setc(c, value); + Sc.setc(c, value); } - return result; - } - gl_spharmonics operator-(gl_spharmonics rhs) { - return (*this) + (rhs * (T)(-1)); - } + 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>& 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>& vlist, size_t nc) { + Sd.pdf(vlist, nc); + Sc = Sd; + } -}; //end gl_spharmonics + void slices(size_t s) { + N = s; + } + size_t slices() { + return N; + } -}; //end namespace stim + void rendermode(bool displace, bool color, bool mag = true) { + displacement = displace; + colormap = color; + magnitude = mag; + } + }; +} -#endif +#endif \ No newline at end of file -- libgit2 0.21.4