codebase / stimlib

Commit 01707489309c7fd07f3372e43063de936141e4cb

Authored by David Mayerich
1 parent d9547496

implemented a new spherical harmonics renderer

Showing 4 changed files with 711 additions and 421 deletions   Show diff stats
Inline Side-by-side
stim/image/image.h
  Show/Hide comments   View file @0170748
@@ -137,10 +137,6 @@ public: @@ -137,10 +137,6 @@ public:
137 std::cout << "Error opening input file in image::load_netpbm()" << std::endl; 137 std::cout << "Error opening input file in image::load_netpbm()" << std::endl;
138 exit(1); 138 exit(1);
139 } 139 }
140 - if (sizeof(T) != 1) {  
141 - std::cout << "Error in image::load_netpbm() - data type must be 8-bit integer." << std::endl;  
142 - exit(1);  
143 - }  
144 140
145 size_t nc; //allocate space for the number of channels 141 size_t nc; //allocate space for the number of channels
146 char format[2]; //allocate space to hold the image format tag 142 char format[2]; //allocate space to hold the image format tag
@@ -187,9 +183,12 @@ public: @@ -187,9 +183,12 @@ public:
187 } 183 }
188 size_t h = atoi(sh.c_str()); //convert the string into an integer 184 size_t h = atoi(sh.c_str()); //convert the string into an integer
189 185
190 - allocate(w, h, nc); //allocate space for the image  
191 - infile.read((char*)img, size()); //copy the binary data from the file to the image  
192 - infile.close(); 186 + allocate(w, h, nc); //allocate space for the image
  187 + unsigned char* buffer = (unsigned char*)malloc(w * h * nc); //create a buffer to store the read data
  188 + infile.read((char*)buffer, size()); //copy the binary data from the file to the image
  189 + infile.close(); //close the file
  190 + for (size_t n = 0; n < size(); n++) img[n] = (T)buffer[n]; //copy the buffer data into the image
  191 + free(buffer); //free the buffer array
193 } 192 }
194 193
195 194
stim/math/spharmonics.h
  Show/Hide comments   View file @0170748
@@ -2,280 +2,348 @@ @@ -2,280 +2,348 @@
2 #define STIM_SPH_HARMONICS 2 #define STIM_SPH_HARMONICS
3 3
4 #include <complex> 4 #include <complex>
5 -#include <stim/math/vector.h>  
6 #include <boost/math/special_functions/spherical_harmonic.hpp> 5 #include <boost/math/special_functions/spherical_harmonic.hpp>
7 #include <stim/math/constants.h> 6 #include <stim/math/constants.h>
8 #include <stim/math/random.h> 7 #include <stim/math/random.h>
9 #include <vector> 8 #include <vector>
10 9
11 #define WIRE_SCALE 1.001 10 #define WIRE_SCALE 1.001
12 -namespace stim{ 11 +namespace stim {
13 12
14 -template<class T>  
15 -class spharmonics{ 13 + template<class T>
  14 + class spharmonics {
16 15
17 -public:  
18 - std::vector<T> C; //list of SH coefficients 16 + public:
  17 + std::vector<T> C; //list of SH coefficients
19 18
20 -protected: 19 + protected:
  20 + unsigned int mcN; //number of Monte-Carlo samples
  21 + unsigned int coeff_1d(unsigned int l, int m) { //convert (l,m) to i (1D coefficient index)
  22 + return pow(l + 1, 2) - (l - m) - 1;
  23 + }
  24 + void coeff_2d(size_t c, unsigned int& l, int& m) { //convert a 1D coefficient index into (l, m)
  25 + l = (unsigned int)ceil(sqrt((double)c + 1)) - 1; //the major index is equal to sqrt(c) - 1
  26 + m = (int)(c - (size_t)(l * l)) - (int)l; //the minor index is calculated by finding the difference
  27 + }
21 28
  29 + public:
  30 + spharmonics() {
  31 + mcN = 0;
  32 + }
  33 + spharmonics(size_t c) : spharmonics() {
  34 + resize(c);
  35 + }
22 36
23 - unsigned int mcN; //number of Monte-Carlo samples 37 + void push(T c) {
  38 + C.push_back(c);
  39 + }
24 40
25 - //calculate the value of the SH basis function (l, m) at (theta, phi)  
26 - //here, theta = [0, PI], phi = [0, 2*PI]  
27 - double SH(int l, int m, double theta, double phi){  
28 - //std::complex<T> result = boost::math::spherical_harmonic(l, m, phi, theta);  
29 - //return result.imag() + result.real();  
30 -  
31 - //this calculation is based on calculating the real spherical harmonics:  
32 - // https://en.wikipedia.org/wiki/Spherical_harmonics#Addition_theorem  
33 - if (m < 0) {  
34 - return sqrt(2.0) * pow(-1, m) * boost::math::spherical_harmonic(l, abs(m), phi, theta).imag(); 41 + void resize(unsigned int n) {
  42 + C.resize(n);
35 } 43 }
36 - else if (m == 0) {  
37 - return boost::math::spherical_harmonic(l, m, phi, theta).real(); 44 +
  45 + void setc(unsigned int l, int m, T value) {
  46 + unsigned int c = coeff_1d(l, m);
  47 + C[c] = value;
38 } 48 }
39 - else {  
40 - return sqrt(2.0) * pow(-1, m) * boost::math::spherical_harmonic(l, m, phi, theta).real(); 49 +
  50 + T getc(unsigned int l, int m) {
  51 + unsigned int c = coeff_1d(l, m);
  52 + return C[c];
41 } 53 }
42 - }  
43 54
44 - unsigned int coeff_1d(unsigned int l, int m){  
45 - return pow(l + 1, 2) - (l - m) - 1;  
46 - } 55 + void setc(unsigned int c, T value) {
  56 + C[c] = value;
  57 + }
47 58
48 - 59 + unsigned int getSize() const {
  60 + return C.size();
  61 + }
49 62
  63 + std::vector<T> getC() const {
  64 + return C;
  65 + }
  66 + //calculate the value of the SH basis function (l, m) at (theta, phi)
  67 + //here, theta = [0, PI], phi = [0, 2*PI]
  68 + T SH(unsigned int l, int m, T theta, T phi) {
  69 + //std::complex<T> result = boost::math::spherical_harmonic(l, m, phi, theta);
  70 + //return result.imag() + result.real();
  71 +
  72 + //this calculation is based on calculating the real spherical harmonics:
  73 + // https://en.wikipedia.org/wiki/Spherical_harmonics#Addition_theorem
  74 + if (m < 0) {
  75 + return sqrt(2.0) * pow(-1, m) * boost::math::spherical_harmonic(l, abs(m), phi, theta).imag();
  76 + }
  77 + else if (m == 0) {
  78 + return boost::math::spherical_harmonic(l, m, phi, theta).real();
  79 + }
  80 + else {
  81 + return sqrt(2.0) * pow(-1, m) * boost::math::spherical_harmonic(l, m, phi, theta).real();
  82 + }
  83 + }
50 84
51 -public:  
52 - spharmonics() {  
53 - mcN = 0;  
54 - }  
55 - spharmonics(size_t c) : spharmonics() {  
56 - resize(c);  
57 - } 85 + /// Calculate the spherical harmonic result given a 1D coefficient index
  86 + T SH(size_t c, T theta, T phi) {
  87 + unsigned int l;
  88 + int m;
  89 + coeff_2d(c, l, m);
  90 + return SH(l, m, theta, phi);
  91 + }
58 92
59 - void push(double c){  
60 - C.push_back(c);  
61 - }  
62 93
63 - void resize(unsigned int n){  
64 - C.resize(n);  
65 - }  
66 94
67 - void setc(unsigned int l, int m, T value){  
68 - unsigned int c = coeff_1d(l, m);  
69 - C[c] = value;  
70 - } 95 + /// Initialize Monte-Carlo sampling of a function using N spherical harmonics coefficients
71 96
72 - void setc(unsigned int c, T value){  
73 - C[c] = value;  
74 - } 97 + /// @param N is the number of spherical harmonics coefficients used to represent the user function
  98 + void mcBegin(unsigned int coefficients) {
  99 + C.resize(coefficients, 0);
  100 + mcN = 0;
  101 + }
75 102
76 - unsigned int getSize() const{  
77 - return C.size();  
78 - } 103 + void mcBegin(unsigned int l, int m) {
  104 + unsigned int c = pow(l + 1, 2) - (l - m);
  105 + mcBegin(c);
  106 + }
79 107
80 - std::vector<T> getC() const{  
81 - return C;  
82 - } 108 + void mcSample(T theta, T phi, T val) {
83 109
84 - 110 + int l, m;
  111 + T sh;
85 112
86 - /// Initialize Monte-Carlo sampling of a function using N spherical harmonics coefficients 113 + l = m = 0;
  114 + for (unsigned int i = 0; i < C.size(); i++) {
87 115
88 - /// @param N is the number of spherical harmonics coefficients used to represent the user function  
89 - void mcBegin(unsigned int coefficients){  
90 - C.resize(coefficients, 0);  
91 - mcN = 0;  
92 - } 116 + sh = SH(l, m, theta, phi);
  117 + C[i] += sh * val;
93 118
94 - void mcBegin(unsigned int l, int m){  
95 - unsigned int c = pow(l + 1, 2) - (l - m);  
96 - mcBegin(c);  
97 - } 119 + m++; //increment m
98 120
99 - void mcSample(double theta, double phi, double val){ 121 + //if we're in a new tier, increment l and set m = -l
  122 + if (m > l) {
  123 + l++;
  124 + m = -l;
  125 + }
  126 + } //end for all coefficients
100 127
101 - int l, m;  
102 - double sh; 128 + //increment the number of samples
  129 + mcN++;
103 130
104 - l = m = 0;  
105 - for(unsigned int i = 0; i < C.size(); i++){ 131 + } //end mcSample()
106 132
107 - sh = SH(l, m, theta, phi);  
108 - C[i] += sh * val; 133 + void mcEnd() {
109 134
110 - m++; //increment m 135 + //divide all coefficients by the number of samples
  136 + for (unsigned int i = 0; i < C.size(); i++)
  137 + C[i] /= mcN;
  138 + }
111 139
112 - //if we're in a new tier, increment l and set m = -l  
113 - if(m > l){  
114 - l++;  
115 - m = -l; 140 + /// Generates a PDF describing the probability distribution of points on a spherical surface
  141 + /// @param sph_pts is a list of points in spherical coordinates (theta, phi) where theta = [0, 2pi] and phi = [0, pi]
  142 + /// @param l is the maximum degree of the spherical harmonic function
  143 + /// @param m is the maximum order
  144 + void pdf(std::vector<stim::vec3<T> > sph_pts, unsigned int l, int m) {
  145 + mcBegin(l, m); //begin spherical harmonic sampling
  146 + unsigned int nP = sph_pts.size();
  147 + for (unsigned int p = 0; p < nP; p++) {
  148 + mcSample(sph_pts[p][1], sph_pts[p][2], 1.0);
116 } 149 }
117 - } //end for all coefficients 150 + mcEnd();
  151 + }
  152 +
  153 + void pdf(std::vector<stim::vec3<T> > sph_pts, size_t c) {
  154 + unsigned int l;
  155 + int m;
  156 + coeff_2d(c, l, m);
  157 + pdf(sph_pts, l, m);
  158 + }
  159 +
  160 + /// Project a set of samples onto a spherical harmonic basis
  161 + void project(std::vector<stim::vec3<T> > sph_pts, unsigned int l, int m) {
  162 + mcBegin(l, m); //begin spherical harmonic sampling
  163 + unsigned int nP = sph_pts.size();
  164 + for (unsigned int p = 0; p < nP; p++) {
  165 + mcSample(sph_pts[p][1], sph_pts[p][2], sph_pts[p][0]);
  166 + }
  167 + mcEnd();
  168 + }
  169 + void project(std::vector<stim::vec3<T> > sph_pts, size_t c) {
  170 + unsigned int l;
  171 + int m;
  172 + coeff_2d(c, l, m);
  173 + project(sph_pts, l, m);
  174 + }
118 175
119 - //increment the number of samples  
120 - mcN++; 176 + /// Generates a PDF describing the density distribution of points on a sphere
  177 + /// @param sph_pts is a list of points in cartesian coordinates
  178 + /// @param l is the maximum degree of the spherical harmonic function
  179 + /// @param m is the maximum order
  180 + /// @param c is the centroid of the points in sph_pts. DEFAULT 0,0,0
  181 + /// @param n is the number of points of the surface of the sphere used to create the PDF. DEFAULT 1000
  182 + /// @param norm, a boolean that sets where the output vectors will be normalized between 0 and 1.
  183 + /// @param
  184 + /*void pdf(std::vector<stim::vec3<T> > sph_pts, unsigned int l, int m, stim::vec3<T> c = stim::vec3<T>(0, 0, 0), unsigned int n = 1000, bool norm = true, std::vector<T> w = std::vector<T>())
  185 + {
  186 + std::vector<double> weights; ///the weight at each point on the surface of the sphere.
  187 + // weights.resize(n);
  188 + unsigned int nP = sph_pts.size();
  189 + std::vector<stim::vec3<T> > sphere = stim::Random<T>::sample_sphere(n, 1.0, stim::TAU);
  190 + if (w.size() < nP)
  191 + w = std::vector<T>(nP, 1.0);
  192 +
  193 + for (int i = 0; i < n; i++)
  194 + {
  195 + T val = 0;
  196 + for (int j = 0; j < nP; j++)
  197 + {
  198 + stim::vec3<T> temp = sph_pts[j] - c;
  199 + if (temp.dot(sphere[i]) > 0)
  200 + val += pow(temp.dot(sphere[i]), 4)*w[j];
  201 + }
  202 + weights.push_back(val);
  203 + }
121 204
122 - } //end mcSample() 205 + mcBegin(l, m); //begin spherical harmonic sampling
123 206
124 - void mcEnd(){ 207 + if (norm)
  208 + {
  209 + T min = *std::min_element(weights.begin(), weights.end());
  210 + T max = *std::max_element(weights.begin(), weights.end());
  211 + for (unsigned int i = 0; i < n; i++)
  212 + {
  213 + stim::vec3<T> sph = sphere[i].cart2sph();
  214 + mcSample(sph[1], sph[2], (weights[i] - min) / (max - min));
  215 + }
125 216
126 - //divide all coefficients by the number of samples  
127 - for(unsigned int i = 0; i < C.size(); i++)  
128 - C[i] /= mcN;  
129 - } 217 + }
  218 + else {
  219 + for (unsigned int i = 0; i < n; i++)
  220 + {
  221 + stim::vec3<T> sph = sphere[i].cart2sph();
  222 + mcSample(sph[1], sph[2], weights[i]);
  223 + }
  224 + }
  225 + mcEnd();
  226 + }*/
130 227
131 - /// Generates a PDF describing the probability distribution of points on a spherical surface  
132 - /// @param sph_pts is a list of points in spherical coordinates (theta, phi) where theta = [0, 2pi] and phi = [0, pi]  
133 - /// @param l is the maximum degree of the spherical harmonic function  
134 - /// @param m is the maximum order 228 + std::string str() {
135 229
136 - void pdf(std::vector<stim::vec<double> > sph_pts, unsigned int l, int m){  
137 -  
138 - mcBegin( l, m ); //begin spherical harmonic sampling 230 + std::stringstream ss;
139 231
140 - unsigned int nP = sph_pts.size(); 232 + int l, m;
  233 + l = m = 0;
  234 + for (unsigned int i = 0; i < C.size(); i++) {
141 235
142 - for(unsigned int p = 0; p < nP; p++){  
143 - mcSample(sph_pts[p][1], sph_pts[p][2], 1.0);  
144 - } 236 + ss << C[i] << '\t';
145 237
146 - mcEnd();  
147 - }  
148 -  
149 - /// Generates a PDF describing the probability distribution of points on a spherical surface  
150 - /// @param sph_pts is a list of points in cartesian coordinates  
151 - /// @param l is the maximum degree of the spherical harmonic function  
152 - /// @param m is the maximum order  
153 - /// @param c is the centroid of the points in sph_pts. DEFAULT 0,0,0  
154 - /// @param n is the number of points of the surface of the sphere used to create the PDF. DEFAULT 1000  
155 - void pdf(std::vector<stim::vec3<double> > sph_pts, unsigned int l, int m, stim::vec3<double> c = stim::vec3<double>(0,0,0), unsigned int n = 1000)  
156 - {  
157 - std::vector<double> weights; ///the weight at each point on the surface of the sphere.  
158 -// weights.resize(n);  
159 - unsigned int nP = sph_pts.size();  
160 - std::vector<stim::vec3<double> > sphere = stim::Random<double>::sample_sphere(n, 1.0, stim::TAU);  
161 - for(int i = 0; i < n; i++)  
162 - {  
163 - double val = 0;  
164 - for(int j = 0; j < nP; j++)  
165 - {  
166 - stim::vec3<double> temp = sph_pts[j] - c;  
167 - if(temp.dot(sphere[i]) > 0)  
168 - val += pow(temp.dot(sphere[i]),4);  
169 - }  
170 - weights.push_back(val);  
171 - }  
172 -  
173 -  
174 - mcBegin(l, m); //begin spherical harmonic sampling  
175 -  
176 - for(unsigned int i = 0; i < n; i++)  
177 - {  
178 - stim::vec3<double> sph = sphere[i].cart2sph();  
179 - mcSample(sph[1], sph[2], weights[i]);  
180 - } 238 + m++; //increment m
181 239
182 - mcEnd();  
183 - } 240 + //if we're in a new tier, increment l and set m = -l
  241 + if (m > l) {
  242 + l++;
  243 + m = -l;
184 244
185 - std::string str(){ 245 + ss << std::endl;
186 246
187 - std::stringstream ss; 247 + }
  248 + }
188 249
189 - int l, m;  
190 - l = m = 0;  
191 - for(unsigned int i = 0; i < C.size(); i++){  
192 -  
193 - ss<<C[i]<<'\t'; 250 + return ss.str();
194 251
195 - m++; //increment m  
196 252
197 - //if we're in a new tier, increment l and set m = -l  
198 - if(m > l){  
199 - l++;  
200 - m = -l; 253 + }
201 254
202 - ss<<std::endl;  
203 - 255 + /// Returns the value of the function at coordinate (theta, phi)
  256 + T p(T theta, T phi) {
  257 + T fx = 0;
  258 +
  259 + int l = 0;
  260 + int m = 0;
  261 + for (unsigned int i = 0; i < C.size(); i++) {
  262 + fx += C[i] * SH(l, m, theta, phi);
  263 + m++;
  264 + if (m > l) {
  265 + l++;
  266 + m = -l;
  267 + }
204 } 268 }
  269 + return fx;
205 } 270 }
206 271
207 - return ss.str(); 272 + /// Returns the derivative of the spherical function with respect to theta
  273 + /// return value is in cartesian coordinates
  274 + vec3<T> dtheta(T theta, T phi, T d = 0.01) {
  275 + T r = p(theta, phi); //calculate the value of the spherical function at three points
  276 + T rt = p(theta + d, phi);
  277 + //double rp = p(theta, phi + d);
208 278
  279 + vec3<T> s(r, theta, phi); //get the spherical coordinate position for all three points
  280 + vec3<T> st(rt, theta + d, phi);
  281 + //vec3<double> sp(rp, theta, phi + d);
209 282
210 - } 283 + vec3<T> c = s.sph2cart();
  284 + vec3<T> ct = st.sph2cart();
  285 + //vec3<double> cp = sp.sph2cart();
211 286
212 - /// Returns the value of the function at the coordinate (theta, phi) 287 + vec3<T> dt = (ct - c)/d; //calculate the derivative
  288 + return dt;
  289 + }
213 290
214 - /// @param theta = [0, 2pi]  
215 - /// @param phi = [0, pi]  
216 - double operator()(double theta, double phi){ 291 + /// Returns the derivative of the spherical function with respect to phi
  292 + /// return value is in cartesian coordinates
  293 + vec3<T> dphi(T theta, T phi, T d = 0.01) {
  294 + T r = p(theta, phi); //calculate the value of the spherical function at three points
  295 + //double rt = p(theta + d, phi);
  296 + T rp = p(theta, phi + d);
217 297
218 - double fx = 0; 298 + vec3<T> s(r, theta, phi); //get the spherical coordinate position for all three points
  299 + //vec3<double> st(rt, theta + d, phi);
  300 + vec3<T> sp(rp, theta, phi + d);
219 301
220 - int l = 0;  
221 - int m = 0;  
222 - for(unsigned int i = 0; i < C.size(); i++){  
223 - fx += C[i] * SH(l, m, theta, phi);  
224 - m++;  
225 - if(m > l){  
226 - l++;  
227 - m = -l;  
228 - } 302 + vec3<T> c = s.sph2cart();
  303 + //vec3<double> ct = st.sph2cart();
  304 + vec3<T> cp = sp.sph2cart();
229 305
  306 + vec3<T> dp = (cp - c) / d; //calculate the derivative
  307 + return dp;
230 } 308 }
231 -  
232 - return fx;  
233 - }  
234 -  
235 - /// Fill an NxN grid with the spherical function for theta = [0 2pi] and phi = [0 pi]  
236 - void grid(T* data, size_t N){  
237 - double dt = stim::TAU / (double)N; //calculate the step size in each direction  
238 - double dp = stim::PI / (double)(N - 1);  
239 - for(size_t ti = 0; ti < N; ti++){  
240 - for(size_t pi = 0; pi < N){  
241 - data[pi * N + ti] = (*this)((double)ti * dt, (double)pi * dp);  
242 - } 309 +
  310 + /// Returns the value of the function at the coordinate (theta, phi)
  311 + /// @param theta = [0, 2pi]
  312 + /// @param phi = [0, pi]
  313 + T operator()(T theta, T phi) {
  314 + return p(theta, phi);
243 } 315 }
244 - }  
245 -/*  
246 - //overload arithmetic operations  
247 - spharmonics<T> operator*(T rhs) const {  
248 - spharmonics<T> result(C.size()); //create a new spherical harmonics object  
249 - for (size_t c = 0; c < C.size(); c++) //for each coefficient  
250 - result.C[c] = C[c] * rhs; //calculate the factor and store the result in the new spharmonics object  
251 - return result;  
252 - }  
253 -  
254 - spharmonics<T> operator+(spharmonics<T> rhs) {  
255 - size_t low = std::min(C.size(), rhs.C.size()); //store the number of coefficients in the lowest object  
256 - size_t high = std::max(C.size(), rhs.C.size()); //store the number of coefficients in the result  
257 - bool rhs_lowest = false; //true if rhs has the lowest number of coefficients  
258 - if (rhs.C.size() < C.size()) rhs_lowest = true; //if rhs has a lower number of coefficients, set the flag  
259 -  
260 - spharmonics<T> result(high); //create a new object  
261 - size_t c;  
262 - for (c = 0; c < low; c++) //perform the first batch of additions  
263 - result.C[c] = C[c] + rhs.C[c]; //perform the addition  
264 -  
265 - for (c = low; c < high; c++) {  
266 - if (rhs_lowest)  
267 - result.C[c] = C[c];  
268 - else  
269 - result.C[c] = rhs.C[c]; 316 + /// Fill an NxN grid with the spherical function for theta = [0 2pi] and phi = [0 pi]
  317 + void get_func(T* data, size_t X, size_t Y) {
  318 + T dt = stim::TAU / (T)X; //calculate the step size in each direction
  319 + T dp = stim::PI / (T)(Y - 1);
  320 + for (size_t ti = 0; ti < X; ti++) {
  321 + for (size_t pi = 0; pi < Y; pi++) {
  322 + data[pi * X + ti] = (*this)((T)ti * dt, (T)pi * dp);
  323 + }
  324 + }
270 } 325 }
271 - return result;  
272 - }  
273 326
274 - spharmonics<T> operator-(spharmonics<T> rhs) {  
275 - return (*this) + (rhs * (T)(-1));  
276 - }  
277 -*/  
278 -}; //end class sph_harmonics 327 + /// Project a spherical function onto the basis using C coefficients
  328 + /// @param data is a pointer to the function values in (theta, phi) coordinates
  329 + /// @param N is the number of samples along each axis, where theta = [0 2pi), phi = [0 pi]
  330 + void project(T* data, size_t x, size_t y, size_t nc) {
  331 + stim::cpu2image(data, "test.ppm", x, y, stim::cmBrewer);
  332 + C.resize(nc, 0); //resize the coefficient array to store the necessary coefficients
  333 + T dtheta = stim::TAU / (T)(x - 1); //calculate the grid spacing along theta
  334 + T dphi = stim::PI / (T)y; //calculate the grid spacing along phi
  335 + T theta, phi;
  336 + for (size_t c = 0; c < nc; c++) { //for each coefficient
  337 + for (size_t theta_i = 0; theta_i < x; theta_i++) { //for each coordinate in the provided array
  338 + theta = theta_i * dtheta; //calculate theta
  339 + for (size_t phi_i = 0; phi_i < y; phi_i++) {
  340 + phi = phi_i * dphi; //calculate phi
  341 + C[c] += data[phi_i * x + theta_i] * SH(c, theta, phi) * dtheta * dphi * sin(phi);
  342 + }
  343 + }
  344 + }
  345 + }
  346 + }; //end class sph_harmonics
279 347
280 348
281 349
@@ -283,4 +351,4 @@ public: @@ -283,4 +351,4 @@ public:
283 } 351 }
284 352
285 353
286 -#endif 354 +#endif
287 \ No newline at end of file 355 \ No newline at end of file
stim/visualization/gl_spharmonics-dep.h 0 → 100644
  Show/Hide comments   View file @0170748
  1 +#ifndef STIM_GL_SPHARMONICS_H
  2 +#define STIM_GL_SPHARMONICS_H
  3 +
  4 +#include <GL/gl.h>
  5 +
  6 +#include <stim/gl/error.h>
  7 +#include <stim/visualization/colormap.h>
  8 +#include <stim/math/spharmonics.h>
  9 +
  10 +namespace stim{
  11 +
  12 +
  13 +template <typename T>
  14 +class gl_spharmonics : public spharmonics<T>{
  15 +
  16 +protected:
  17 +
  18 + using stim::spharmonics<T>::SH;
  19 + stim::spharmonics<T> surface; //harmonic that stores the surface information
  20 + stim::spharmonics<T> color; //harmonic that stores the color information
  21 + T* func_surface; //stores the raw function data (samples at each point)
  22 + T* func_color; //stores the raw color data (samples at each point)
  23 + GLuint color_tex; //texture map that acts as a colormap for the spherical function
  24 + unsigned int N; //resolution of the spherical grid
  25 +
  26 + ///generates the
  27 + void gen_surface(){
  28 + //allocate memory and initialize the function to zero
  29 + func_surface = (T*) malloc(N * N * sizeof(T));
  30 + memset(func_surface, 0, sizeof(T) * N * N);
  31 +
  32 + double theta, phi;
  33 + double result;
  34 + int l, m;
  35 +
  36 +
  37 + l = m = 0;
  38 + //for each coefficient
  39 + for(unsigned int c = 0; c < surface.C.size(); c++){
  40 +
  41 + //iterate through the entire 2D grid representing the function
  42 + for(unsigned int xi = 0; xi < N; xi++){
  43 + for(unsigned int yi = 0; yi < N; yi++){
  44 +
  45 + //get the spherical coordinates for each grid point
  46 + theta = (2 * PI) * ((double)xi / (N-1));
  47 + phi = PI * ((double)yi / (N-1));
  48 +
  49 + //sum the contribution of the current spherical harmonic based on the coefficient
  50 + result = surface.C[c] * SH(l, m, theta, phi);
  51 +
  52 + //store the result in a 2D array (which will later be used as a texture map)
  53 + func_surface[yi * N + xi] += result;
  54 + }
  55 + }
  56 +
  57 + //keep track of m and l here
  58 + m++; //increment m
  59 + //if we're in a new tier, increment l and set m = -l
  60 + if(m > l){
  61 + l++;
  62 + m = -l;
  63 + }
  64 + }
  65 + }
  66 +
  67 + ///generates the color function
  68 + void gen_color(){
  69 +
  70 + //initialize the function to zero
  71 + func_color = (T*) malloc(N * N * sizeof(T));
  72 + memset(func_color, 0, sizeof(T) * N * N);
  73 +
  74 + double theta, phi;
  75 + double result;
  76 + int l, m;
  77 +
  78 +
  79 + l = m = 0;
  80 + //for each coefficient
  81 + for(unsigned int c = 0; c < color.C.size(); c++){
  82 +
  83 + //iterate through the entire 2D grid representing the function
  84 + for(unsigned int xi = 0; xi < N; xi++){
  85 + for(unsigned int yi = 0; yi < N; yi++){
  86 +
  87 + //get the spherical coordinates for each grid point
  88 + theta = (2 * PI) * ((double)xi / (N-1));
  89 + phi = PI * ((double)yi / (N-1));
  90 +
  91 + //sum the contribution of the current spherical harmonic based on the coefficient
  92 + result = color.C[c] * SH(l, m, theta, phi);
  93 +
  94 + //store the result in a 2D array (which will later be used as a texture map)
  95 + func_color[yi * N + xi] += result;
  96 + }
  97 + }
  98 +
  99 + //keep track of m and l here
  100 + m++; //increment m
  101 + //if we're in a new tier, increment l and set m = -l
  102 + if(m > l){
  103 + l++;
  104 + m = -l;
  105 + }
  106 + }
  107 + }
  108 +
  109 + void gl_prep_draw(){
  110 +
  111 + //enable depth testing
  112 + //this has to be used instead of culling because the sphere can have negative values
  113 + glEnable(GL_DEPTH_TEST);
  114 + glDepthMask(GL_TRUE);
  115 + glEnable(GL_TEXTURE_2D); //enable 2D texture mapping
  116 + }
  117 +
  118 + //draw a texture mapped sphere representing the function surface
  119 + void gl_draw_sphere() {
  120 +
  121 + //bind the 2D texture representing the color map
  122 + glBindTexture(GL_TEXTURE_2D, color_tex);
  123 +
  124 + //Draw the Sphere
  125 + int i, j;
  126 +
  127 + for(i = 1; i <= N-1; i++) {
  128 + double phi0 = PI * ((double) (i - 1) / (N-1));
  129 + double phi1 = PI * ((double) i / (N-1));
  130 +
  131 + glBegin(GL_QUAD_STRIP);
  132 + for(j = 0; j <= N; j++) {
  133 +
  134 + //calculate the indices into the function array
  135 + int phi0_i = i-1;
  136 + int phi1_i = i;
  137 + int theta_i = j;
  138 + if(theta_i == N)
  139 + theta_i = 0;
  140 +
  141 + double v0 = func_surface[phi0_i * N + theta_i];
  142 + double v1 = func_surface[phi1_i * N + theta_i];
  143 +
  144 + v0 = fabs(v0);
  145 + v1 = fabs(v1);
  146 +
  147 +
  148 + double theta = 2 * PI * (double) (j - 1) / N;
  149 + double x0 = v0 * cos(theta) * sin(phi0);
  150 + double y0 = v0 * sin(theta) * sin(phi0);
  151 + double z0 = v0 * cos(phi0);
  152 +
  153 + double x1 = v1 * cos(theta) * sin(phi1);
  154 + double y1 = v1 * sin(theta) * sin(phi1);
  155 + double z1 = v1 * cos(phi1);
  156 +
  157 + glTexCoord2f(theta / (2 * PI), phi0 / PI);
  158 + glVertex3f(x0, y0, z0);
  159 +
  160 + glTexCoord2f(theta / (2 * PI), phi1 / PI);
  161 + glVertex3f(x1, y1, z1);
  162 + }
  163 + glEnd();
  164 + }
  165 +
  166 + glBindTexture(GL_TEXTURE_2D, 0);
  167 + }
  168 +
  169 + void gen_texture()
  170 + {
  171 + //allocate space for the color map
  172 + unsigned int bytes = N * N * sizeof(unsigned char) * 3;
  173 + unsigned char* color_image;
  174 + color_image = (unsigned char*) malloc(bytes);
  175 +
  176 + //generate a colormap from the function
  177 + stim::cpu2cpu<double>(func_color, color_image, N*N, stim::cmBrewer);
  178 +
  179 + //prep everything for drawing
  180 + gl_prep_draw();
  181 +
  182 + //generate an OpenGL texture map in the current context
  183 + glGenTextures(1, &color_tex);
  184 + //bind the texture
  185 + glBindTexture(GL_TEXTURE_2D, color_tex);
  186 +
  187 + //copy the color data from the buffer to the GPU
  188 + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, color_image);
  189 +
  190 + //initialize all of the texture parameters
  191 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
  192 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
  193 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
  194 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  195 + glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
  196 +
  197 + //free the buffer
  198 + free(color_image);
  199 + }
  200 +
  201 + void gl_init(unsigned int n){
  202 +
  203 + //set the sphere resolution
  204 + N = n;
  205 +
  206 + //set up the surface data.
  207 + gen_surface();
  208 +
  209 + //set up the color data and generate the corresponding texture.
  210 + gen_color();
  211 + gen_texture();
  212 +
  213 + }
  214 +
  215 +public:
  216 + gl_spharmonics<T>(unsigned int n = 256) {
  217 + gl_init(n);
  218 + }
  219 +
  220 + gl_spharmonics<T>( stim::spharmonics<T> in_function, unsigned int n = 256)
  221 + {
  222 + surface = in_function;
  223 + color = in_function;
  224 + gl_init(n);
  225 + }
  226 +
  227 + gl_spharmonics<T>( stim::spharmonics<T> in_function, stim::spharmonics<T> in_color, unsigned int n = 256)
  228 + {
  229 + surface = in_function;
  230 + color = in_color;
  231 + gl_init(n);
  232 + }
  233 +
  234 + gl_spharmonics<T>& operator=(const spharmonics<T> rhs) {
  235 + gl_spharmonics<T> result(rhs.C.size());
  236 + result.C = spharmonics<T>::rhs.C;
  237 + return result;
  238 + }
  239 +
  240 + void glRender(){
  241 + //set all OpenGL parameters required for drawing
  242 + gl_prep_draw();
  243 +
  244 + //draw the sphere
  245 + gl_draw_sphere();
  246 + }
  247 +
  248 + void glInit(unsigned int n = 256){
  249 + gl_init(n);
  250 + }
  251 +}; //end gl_spharmonics
  252 +
  253 +
  254 +}; //end namespace stim
  255 +
  256 +
  257 +
  258 +
  259 +#endif
stim/visualization/gl_spharmonics.h
  Show/Hide comments   View file @0170748
1 #ifndef STIM_GL_SPHARMONICS_H 1 #ifndef STIM_GL_SPHARMONICS_H
2 #define STIM_GL_SPHARMONICS_H 2 #define STIM_GL_SPHARMONICS_H
3 3
4 -#include <GL/gl.h>  
5 - 4 +#include <stim/math/spharmonics.h>
6 #include <stim/gl/error.h> 5 #include <stim/gl/error.h>
  6 +#include <stim/math/vec3.h>
  7 +#include <stim/math/constants.h>
7 #include <stim/visualization/colormap.h> 8 #include <stim/visualization/colormap.h>
8 -#include <stim/math/spharmonics.h>  
9 -  
10 -namespace stim{  
11 -  
12 -  
13 -template <typename T>  
14 -class gl_spharmonics : public spharmonics<T>{  
15 -  
16 -protected:  
17 - using spharmonics<T>::C;  
18 - using spharmonics<T>::SH;  
19 -  
20 - T* func; //stores the raw function data (samples at each point)  
21 -  
22 - GLuint color_tex; //texture map that acts as a colormap for the spherical function  
23 -  
24 - unsigned int N; //resolution of the spherical grid  
25 -  
26 - void gen_function(){  
27 -  
28 - //initialize the function to zero  
29 - memset(func, 0, sizeof(double) * N * N);  
30 -  
31 - double theta, phi;  
32 - double result;  
33 - int l, m;  
34 9
  10 +namespace stim {
  11 +
  12 + template<typename T>
  13 + class gl_spharmonics {
  14 +
  15 + GLuint dlist;
  16 + GLuint tex;
  17 +
  18 + bool displacement;
  19 + bool colormap;
  20 + bool magnitude;
  21 +
  22 + void init_tex() {
  23 + T* sfunc = (T*)malloc(N * N * sizeof(T)); //create a 2D array to store the spherical function
  24 + Sc.get_func(sfunc, N, N); //generate the spherical function based on the Sc coefficients
  25 + unsigned char* tex_buffer = (unsigned char*)malloc(3 * N * N); //create a buffer to store the texture map
  26 + stim::cpu2cpu<T>(sfunc, tex_buffer, N * N, stim::cmBrewer); //create a Brewer colormap based on the spherical function
  27 + stim::buffer2image(tex_buffer, "sfunc.ppm", N, N);
  28 +
  29 + if (tex) glDeleteTextures(1, &tex); //if a texture already exists, delete it
  30 + glGenTextures(1, &tex); //create a new texture and store the ID
  31 + glBindTexture(GL_TEXTURE_2D, tex); //bind the texture
  32 + 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
  33 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT); //initialize all of the texture parameters
  34 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
  35 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
  36 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  37 + glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
  38 + }
35 39
36 - l = m = 0;  
37 - //for each coefficient  
38 - for(unsigned int c = 0; c < C.size(); c++){ 40 + public:
  41 + stim::spharmonics<T> Sc; //spherical harmonic representing the color component
  42 + stim::spharmonics<T> Sd; //spherical harmonic representing the displacement component
  43 + size_t N;
  44 +
  45 + gl_spharmonics(size_t slices) {
  46 + N = slices;
  47 + dlist = 0; //initialize the display list index to zero (no list)
  48 + tex = 0; //initialize the texture index to zero (no texture)
  49 + displacement = true;
  50 + colormap = true;
  51 + magnitude = true;
  52 + }
39 53
40 - //iterate through the entire 2D grid representing the function  
41 - for(unsigned int xi = 0; xi < N; xi++){  
42 - for(unsigned int yi = 0; yi < N; yi++){ 54 + ~gl_spharmonics() {
  55 + if (dlist) glDeleteLists(dlist, 1); //delete the display list when the object is destroyed
  56 + }
43 57
44 - //get the spherical coordinates for each grid point  
45 - theta = (2 * PI) * ((double)xi / (N-1));  
46 - phi = PI * ((double)yi / (N-1)); 58 + /// This function renders the spherical harmonic to the current OpenGL context
  59 + void render() {
  60 + //glShadeModel(GL_FLAT);
  61 + glPushAttrib(GL_ENABLE_BIT);
  62 + glDisable(GL_CULL_FACE);
  63 +
  64 + if (!tex) {
  65 + init_tex();
  66 + }
47 67
48 - //sum the contribution of the current spherical harmonic based on the coefficient  
49 - result = C[c] * SH(l, m, theta, phi); 68 + if (colormap) {
  69 + glEnable(GL_TEXTURE_2D);
  70 + glBindTexture(GL_TEXTURE_2D, tex);
  71 + }
50 72
51 - //store the result in a 2D array (which will later be used as a texture map)  
52 - func[yi * N + xi] += result; 73 + if (!dlist) {
  74 + dlist = glGenLists(1);
  75 + glNewList(dlist, GL_COMPILE);
  76 + glPushAttrib(GL_ENABLE_BIT);
  77 + glEnable(GL_NORMALIZE);
  78 +
  79 + //Draw the Sphere
  80 + size_t theta_i, phi_i;
  81 + T d_theta = (T)stim::TAU / (T)N;
  82 + T d_phi = (T)stim::PI / (T)(N-1);
  83 +
  84 + for (phi_i = 1; phi_i < N; phi_i++) {
  85 + T phi = phi_i * d_phi;
  86 + glBegin(GL_QUAD_STRIP);
  87 + for (theta_i = 0; theta_i <= N; theta_i++) {
  88 + T theta = (N - theta_i) * d_theta;
  89 + float theta_t = 1 - (float)theta_i / (float)N;
  90 +
  91 + T r;
  92 + if (!displacement) r = 1; //if no displacement, set the r value to 1 (renders a unit sphere)
  93 + else r = Sd(theta, phi); //otherwise calculate the displacement value
  94 +
  95 + glColor3f(1.0f, 1.0f, 1.0f);
  96 + if (!colormap) { //if no colormap is being rendered
  97 + if (r < 0) glColor3f(1.0, 0.0, 0.0); //if r is negative, render it red
  98 + else glColor3f(0.0, 1.0, 0.0); //otherwise render in green
  99 + }
  100 + if (magnitude) { //if the magnitude is being displaced, calculate the magnitude of r
  101 + if (r < 0) r = -r;
  102 + }
  103 + stim::vec3<T> s(r, theta, phi);
  104 + stim::vec3<T> c = s.sph2cart();
  105 + stim::vec3<T> n; //allocate a value to store the normal
  106 + if (!displacement) n = c; //if there is no displacement, the normal is spherical
  107 + else n = Sd.dphi(theta, phi).cross(Sd.dtheta(theta, phi)); //otherwise calculate the normal as the cross product of derivatives
  108 +
  109 + glTexCoord2f(theta_t, (float)phi_i / (float)N);
  110 + //std::cout << theta_t <<" "<<(float)phi_i / (float)N << "----------------";
  111 + glNormal3f(n[0], n[1], n[2]);
  112 + glVertex3f(c[0], c[1], c[2]);
  113 +
  114 + T r1;
  115 + if (!displacement) r1 = 1;
  116 + else r1 = Sd(theta, phi - d_phi);
  117 + if (!colormap) { //if no colormap is being rendered
  118 + if (r1 < 0) glColor3f(1.0, 0.0, 0.0); //if r1 is negative, render it red
  119 + else glColor3f(0.0, 1.0, 0.0); //otherwise render in green
  120 + }
  121 + if (magnitude) { //if the magnitude is being rendered, calculate the magnitude of r
  122 + if (r1 < 0) r1 = -r1;
  123 + }
  124 + stim::vec3<T> s1(r1, theta, phi - d_phi);
  125 + stim::vec3<T> c1 = s1.sph2cart();
  126 + stim::vec3<T> n1;
  127 + if (!displacement) n1 = c1;
  128 + else n1 = Sd.dphi(theta, phi - d_phi).cross(Sd.dtheta(theta, phi - d_phi));
  129 +
  130 + //std::cout << theta_t << " " << (float)(phi_i - 1) / (float)N << std::endl;
  131 + glTexCoord2f(theta_t, 1.0/(2*(N)) + (float)(phi_i-1) / (float)N);
  132 + glNormal3f(n1[0], n1[1], n1[2]);
  133 + glVertex3f(c1[0], c1[1], c1[2]);
  134 + }
  135 + glEnd();
53 } 136 }
  137 + glPopAttrib();
  138 + glEndList();
54 } 139 }
  140 + glCallList(dlist); //call the display list to render
55 141
56 - //keep track of m and l here  
57 - m++; //increment m  
58 - //if we're in a new tier, increment l and set m = -l  
59 - if(m > l){  
60 - l++;  
61 - m = -l;  
62 - } 142 + glPopAttrib();
63 } 143 }
64 - }  
65 -  
66 - void gl_prep_draw(){  
67 -  
68 - //enable depth testing  
69 - //this has to be used instead of culling because the sphere can have negative values  
70 - glEnable(GL_DEPTH_TEST);  
71 - glDepthMask(GL_TRUE);  
72 - glEnable(GL_TEXTURE_2D); //enable 2D texture mapping  
73 - }  
74 -  
75 - //draw a texture mapped sphere representing the function surface  
76 - void gl_draw_sphere() {  
77 -  
78 - //bind the 2D texture representing the color map  
79 - glBindTexture(GL_TEXTURE_2D, color_tex);  
80 -  
81 - //Draw the Sphere  
82 - int i, j;  
83 -  
84 - for(i = 1; i <= N-1; i++) {  
85 - double phi0 = PI * ((double) (i - 1) / (N-1));  
86 - double phi1 = PI * ((double) i / (N-1));  
87 -  
88 - glBegin(GL_QUAD_STRIP);  
89 - for(j = 0; j <= N; j++) {  
90 144
91 - //calculate the indices into the function array  
92 - int phi0_i = i-1;  
93 - int phi1_i = i;  
94 - int theta_i = j;  
95 - if(theta_i == N)  
96 - theta_i = 0;  
97 -  
98 - double v0 = func[phi0_i * N + theta_i];  
99 - double v1 = func[phi1_i * N + theta_i];  
100 -  
101 - v0 = fabs(v0);  
102 - v1 = fabs(v1);  
103 -  
104 -  
105 - double theta = 2 * PI * (double) (j - 1) / N;  
106 - double x0 = v0 * cos(theta) * sin(phi0);  
107 - double y0 = v0 * sin(theta) * sin(phi0);  
108 - double z0 = v0 * cos(phi0);  
109 -  
110 - double x1 = v1 * cos(theta) * sin(phi1);  
111 - double y1 = v1 * sin(theta) * sin(phi1);  
112 - double z1 = v1 * cos(phi1);  
113 -  
114 - glTexCoord2f(theta / (2 * PI), phi0 / PI);  
115 - glVertex3f(x0, y0, z0);  
116 -  
117 - glTexCoord2f(theta / (2 * PI), phi1 / PI);  
118 - glVertex3f(x1, y1, z1);  
119 - }  
120 - glEnd(); 145 + /// Push a coefficient to the spherical harmonic - by default, push applies the component to both the displacement and color SH
  146 + void push(T coeff) {
  147 + Sd.push(coeff);
  148 + Sc.push(coeff);
121 } 149 }
122 - }  
123 -  
124 - void gl_init(unsigned int n){  
125 -  
126 - //set the sphere resolution  
127 - N = n;  
128 -  
129 - //allocate space for the color map  
130 - unsigned int bytes = N * N * sizeof(unsigned char) * 3;  
131 - unsigned char* color_image;  
132 - color_image = (unsigned char*) malloc(bytes);  
133 -  
134 - //allocate space for the function  
135 - func = (double*) malloc(N * N * sizeof(double));  
136 150
137 - //generate a functional representation that will be used for the texture map and vertices  
138 - gen_function();  
139 -  
140 - //generate a colormap from the function  
141 - stim::cpu2cpu<double>(func, color_image, N*N, stim::cmBrewer);  
142 -  
143 - //prep everything for drawing  
144 - gl_prep_draw();  
145 -  
146 - //generate an OpenGL texture map in the current context  
147 - glGenTextures(1, &color_tex);  
148 - //bind the texture  
149 - glBindTexture(GL_TEXTURE_2D, color_tex);  
150 -  
151 - //copy the color data from the buffer to the GPU  
152 - glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, color_image); 151 + /// Resize the spherical harmonic coefficient array
  152 + void resize(size_t s) {
  153 + Sd.resize(s);
  154 + Sc.resize(s);
  155 + }
153 156
154 - //initialize all of the texture parameters  
155 - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);  
156 - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);  
157 - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);  
158 - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);  
159 - glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);  
160 -  
161 - //free the buffer  
162 - free(color_image);  
163 - }  
164 -  
165 -public:  
166 - gl_spharmonics<T>() {  
167 -  
168 - }  
169 -  
170 - gl_spharmonics<T>(const spharmonics<T> copy) : gl_spharmonics<T>() {  
171 - C = copy.C;  
172 - }  
173 -  
174 - gl_spharmonics<T>& operator=(const spharmonics<T> rhs) {  
175 - gl_spharmonics<T> result(rhs.C.size());  
176 - result.C = spharmonics<T>::rhs.C;  
177 - return result;  
178 - }  
179 -  
180 - void glRender(){  
181 - //set all OpenGL parameters required for drawing  
182 - gl_prep_draw();  
183 -  
184 - //draw the sphere  
185 - gl_draw_sphere();  
186 - }  
187 -  
188 - void glInit(unsigned int n = 256){  
189 - gl_init(n);  
190 - gen_function();  
191 - }  
192 -  
193 - //overload arithmetic operations  
194 - gl_spharmonics<T> operator*(T rhs) const {  
195 - gl_spharmonics<T> result(C.size()); //create a new sph erical harmonics object  
196 - for (size_t c = 0; c < C.size(); c++) //for each coefficient  
197 - result.C[c] = C[c] * rhs; //calculat e the factor and store the result in the new spharmonics object  
198 - return result;  
199 - }  
200 -  
201 - gl_spharmonics<T> operator+(gl_spharmonics<T> rhs) {  
202 - size_t low = std::min(C.size(), rhs.C.size()); //store th e number of coefficients in the lowest object  
203 - size_t high = std::max(C.size(), rhs.C.size()); //store th e number of coefficients in the result  
204 - bool rhs_lowest = false; //true if rhs has the lowest number of coefficients  
205 - if (rhs.C.size() < C.size()) rhs_lowest = true; //if rhs has a low er number of coefficients, set the flag  
206 -  
207 - gl_spharmonics<T> result(high); //create a new object  
208 - size_t c;  
209 - for (c = 0; c < low; c++) //perform the first batch of additions  
210 - result.C[c] = C[c] + rhs.C[c]; // perform the addition  
211 -  
212 - for (c = low; c < high; c++) {  
213 - if (rhs_lowest)  
214 - result.C[c] = C[c];  
215 - else  
216 - result.C[c] = rhs.C[c]; 157 + /// Set a spherical harmonic coefficient to the given value
  158 + void setc(unsigned int c, T value) {
  159 + Sd.setc(c, value);
  160 + Sc.setc(c, value);
217 } 161 }
218 - return result;  
219 - }  
220 162
221 - gl_spharmonics<T> operator-(gl_spharmonics<T> rhs) {  
222 - return (*this) + (rhs * (T)(-1));  
223 - } 163 + void project(T* data, size_t x, size_t y, size_t nc) {
  164 + Sd.project(data, x, y, nc);
  165 + Sc = Sd;
  166 + }
224 167
  168 + /// Project a set of samples onto the basis
  169 + void project(std::vector<vec3<float>>& vlist, size_t nc) {
  170 + Sd.project(vlist, nc);
  171 + Sc = Sd;
  172 + }
225 173
  174 + /// Calculate a density function from a list of points in spherical coordinates
  175 + void pdf(std::vector<stim::vec3<T>>& vlist, size_t nc) {
  176 + Sd.pdf(vlist, nc);
  177 + Sc = Sd;
  178 + }
226 179
227 -}; //end gl_spharmonics 180 + void slices(size_t s) {
  181 + N = s;
  182 + }
228 183
  184 + size_t slices() {
  185 + return N;
  186 + }
229 187
230 -}; //end namespace stim 188 + void rendermode(bool displace, bool color, bool mag = true) {
  189 + displacement = displace;
  190 + colormap = color;
  191 + magnitude = mag;
  192 + }
231 193
  194 + };
  195 +}
232 196
233 197
234 198
235 -#endif 199 +#endif
236 \ No newline at end of file 200 \ No newline at end of file