487a9b49
David Mayerich
added the ability...
|
1
2
3
|
#ifndef STIM_SPH_HARMONICS
#define STIM_SPH_HARMONICS
|
9bfae4d0
Pavel Govyadinov
extended gl_sphar...
|
4
|
#include <complex>
|
487a9b49
David Mayerich
added the ability...
|
5
|
#include <boost/math/special_functions/spherical_harmonic.hpp>
|
9bfae4d0
Pavel Govyadinov
extended gl_sphar...
|
6
|
#include <stim/math/constants.h>
|
2e325f29
Pavel Govyadinov
Fixed the pdf imp...
|
7
|
#include <stim/math/random.h>
|
487a9b49
David Mayerich
added the ability...
|
8
9
|
#include <vector>
|
487a9b49
David Mayerich
added the ability...
|
10
|
#define WIRE_SCALE 1.001
|
01707489
David Mayerich
implemented a new...
|
11
|
namespace stim {
|
487a9b49
David Mayerich
added the ability...
|
12
|
|
01707489
David Mayerich
implemented a new...
|
13
14
|
template<class T>
class spharmonics {
|
487a9b49
David Mayerich
added the ability...
|
15
|
|
01707489
David Mayerich
implemented a new...
|
16
17
|
public:
std::vector<T> C; //list of SH coefficients
|
9bfae4d0
Pavel Govyadinov
extended gl_sphar...
|
18
|
|
01707489
David Mayerich
implemented a new...
|
19
20
21
22
23
24
25
26
27
|
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
}
|
487a9b49
David Mayerich
added the ability...
|
28
|
|
01707489
David Mayerich
implemented a new...
|
29
30
31
32
33
34
35
|
public:
spharmonics() {
mcN = 0;
}
spharmonics(size_t c) : spharmonics() {
resize(c);
}
|
487a9b49
David Mayerich
added the ability...
|
36
|
|
01707489
David Mayerich
implemented a new...
|
37
38
39
|
void push(T c) {
C.push_back(c);
}
|
487a9b49
David Mayerich
added the ability...
|
40
|
|
01707489
David Mayerich
implemented a new...
|
41
42
|
void resize(unsigned int n) {
C.resize(n);
|
e3b53161
David Mayerich
fixed real/comple...
|
43
|
}
|
01707489
David Mayerich
implemented a new...
|
44
45
46
47
|
void setc(unsigned int l, int m, T value) {
unsigned int c = coeff_1d(l, m);
C[c] = value;
|
e3b53161
David Mayerich
fixed real/comple...
|
48
|
}
|
01707489
David Mayerich
implemented a new...
|
49
50
51
52
|
T getc(unsigned int l, int m) {
unsigned int c = coeff_1d(l, m);
return C[c];
|
e3b53161
David Mayerich
fixed real/comple...
|
53
|
}
|
487a9b49
David Mayerich
added the ability...
|
54
|
|
01707489
David Mayerich
implemented a new...
|
55
56
57
|
void setc(unsigned int c, T value) {
C[c] = value;
}
|
487a9b49
David Mayerich
added the ability...
|
58
|
|
01707489
David Mayerich
implemented a new...
|
59
60
61
|
unsigned int getSize() const {
return C.size();
}
|
487a9b49
David Mayerich
added the ability...
|
62
|
|
01707489
David Mayerich
implemented a new...
|
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
|
std::vector<T> 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<T> 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();
}
}
|
487a9b49
David Mayerich
added the ability...
|
84
|
|
01707489
David Mayerich
implemented a new...
|
85
86
87
88
89
90
91
|
/// 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);
}
|
487a9b49
David Mayerich
added the ability...
|
92
|
|
487a9b49
David Mayerich
added the ability...
|
93
|
|
487a9b49
David Mayerich
added the ability...
|
94
|
|
01707489
David Mayerich
implemented a new...
|
95
|
/// Initialize Monte-Carlo sampling of a function using N spherical harmonics coefficients
|
487a9b49
David Mayerich
added the ability...
|
96
|
|
01707489
David Mayerich
implemented a new...
|
97
98
99
100
101
|
/// @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;
}
|
487a9b49
David Mayerich
added the ability...
|
102
|
|
01707489
David Mayerich
implemented a new...
|
103
104
105
106
|
void mcBegin(unsigned int l, int m) {
unsigned int c = pow(l + 1, 2) - (l - m);
mcBegin(c);
}
|
9bfae4d0
Pavel Govyadinov
extended gl_sphar...
|
107
|
|
01707489
David Mayerich
implemented a new...
|
108
|
void mcSample(T theta, T phi, T val) {
|
9bfae4d0
Pavel Govyadinov
extended gl_sphar...
|
109
|
|
01707489
David Mayerich
implemented a new...
|
110
111
|
int l, m;
T sh;
|
9bfae4d0
Pavel Govyadinov
extended gl_sphar...
|
112
|
|
01707489
David Mayerich
implemented a new...
|
113
114
|
l = m = 0;
for (unsigned int i = 0; i < C.size(); i++) {
|
487a9b49
David Mayerich
added the ability...
|
115
|
|
01707489
David Mayerich
implemented a new...
|
116
117
|
sh = SH(l, m, theta, phi);
C[i] += sh * val;
|
487a9b49
David Mayerich
added the ability...
|
118
|
|
01707489
David Mayerich
implemented a new...
|
119
|
m++; //increment m
|
487a9b49
David Mayerich
added the ability...
|
120
|
|
01707489
David Mayerich
implemented a new...
|
121
122
123
124
125
126
|
//if we're in a new tier, increment l and set m = -l
if (m > l) {
l++;
m = -l;
}
} //end for all coefficients
|
487a9b49
David Mayerich
added the ability...
|
127
|
|
01707489
David Mayerich
implemented a new...
|
128
129
|
//increment the number of samples
mcN++;
|
487a9b49
David Mayerich
added the ability...
|
130
|
|
01707489
David Mayerich
implemented a new...
|
131
|
} //end mcSample()
|
487a9b49
David Mayerich
added the ability...
|
132
|
|
01707489
David Mayerich
implemented a new...
|
133
|
void mcEnd() {
|
487a9b49
David Mayerich
added the ability...
|
134
|
|
01707489
David Mayerich
implemented a new...
|
135
136
137
138
|
//divide all coefficients by the number of samples
for (unsigned int i = 0; i < C.size(); i++)
C[i] /= mcN;
}
|
487a9b49
David Mayerich
added the ability...
|
139
|
|
01707489
David Mayerich
implemented a new...
|
140
141
142
143
144
145
146
147
148
|
/// 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<stim::vec3<T> > 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);
|
487a9b49
David Mayerich
added the ability...
|
149
|
}
|
01707489
David Mayerich
implemented a new...
|
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
|
mcEnd();
}
void pdf(std::vector<stim::vec3<T> > 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<stim::vec3<T> > 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<stim::vec3<T> > sph_pts, size_t c) {
unsigned int l;
int m;
coeff_2d(c, l, m);
project(sph_pts, l, m);
}
|
487a9b49
David Mayerich
added the ability...
|
175
|
|
01707489
David Mayerich
implemented a new...
|
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
|
/// 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<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>())
{
std::vector<double> weights; ///the weight at each point on the surface of the sphere.
// weights.resize(n);
unsigned int nP = sph_pts.size();
std::vector<stim::vec3<T> > sphere = stim::Random<T>::sample_sphere(n, 1.0, stim::TAU);
if (w.size() < nP)
w = std::vector<T>(nP, 1.0);
for (int i = 0; i < n; i++)
{
T val = 0;
for (int j = 0; j < nP; j++)
{
stim::vec3<T> temp = sph_pts[j] - c;
if (temp.dot(sphere[i]) > 0)
val += pow(temp.dot(sphere[i]), 4)*w[j];
}
weights.push_back(val);
}
|
487a9b49
David Mayerich
added the ability...
|
204
|
|
01707489
David Mayerich
implemented a new...
|
205
|
mcBegin(l, m); //begin spherical harmonic sampling
|
487a9b49
David Mayerich
added the ability...
|
206
|
|
01707489
David Mayerich
implemented a new...
|
207
208
209
210
211
212
213
214
215
|
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<T> sph = sphere[i].cart2sph();
mcSample(sph[1], sph[2], (weights[i] - min) / (max - min));
}
|
487a9b49
David Mayerich
added the ability...
|
216
|
|
01707489
David Mayerich
implemented a new...
|
217
218
219
220
221
222
223
224
225
226
|
}
else {
for (unsigned int i = 0; i < n; i++)
{
stim::vec3<T> sph = sphere[i].cart2sph();
mcSample(sph[1], sph[2], weights[i]);
}
}
mcEnd();
}*/
|
487a9b49
David Mayerich
added the ability...
|
227
|
|
01707489
David Mayerich
implemented a new...
|
228
|
std::string str() {
|
2e325f29
Pavel Govyadinov
Fixed the pdf imp...
|
229
|
|
01707489
David Mayerich
implemented a new...
|
230
|
std::stringstream ss;
|
487a9b49
David Mayerich
added the ability...
|
231
|
|
01707489
David Mayerich
implemented a new...
|
232
233
234
|
int l, m;
l = m = 0;
for (unsigned int i = 0; i < C.size(); i++) {
|
487a9b49
David Mayerich
added the ability...
|
235
|
|
01707489
David Mayerich
implemented a new...
|
236
|
ss << C[i] << '\t';
|
487a9b49
David Mayerich
added the ability...
|
237
|
|
01707489
David Mayerich
implemented a new...
|
238
|
m++; //increment m
|
2e325f29
Pavel Govyadinov
Fixed the pdf imp...
|
239
|
|
01707489
David Mayerich
implemented a new...
|
240
241
242
243
|
//if we're in a new tier, increment l and set m = -l
if (m > l) {
l++;
m = -l;
|
2e325f29
Pavel Govyadinov
Fixed the pdf imp...
|
244
|
|
01707489
David Mayerich
implemented a new...
|
245
|
ss << std::endl;
|
487a9b49
David Mayerich
added the ability...
|
246
|
|
01707489
David Mayerich
implemented a new...
|
247
248
|
}
}
|
487a9b49
David Mayerich
added the ability...
|
249
|
|
01707489
David Mayerich
implemented a new...
|
250
|
return ss.str();
|
487a9b49
David Mayerich
added the ability...
|
251
|
|
487a9b49
David Mayerich
added the ability...
|
252
|
|
01707489
David Mayerich
implemented a new...
|
253
|
}
|
487a9b49
David Mayerich
added the ability...
|
254
|
|
01707489
David Mayerich
implemented a new...
|
255
256
257
258
259
260
261
262
263
264
265
266
267
|
/// Returns the value of the function at coordinate (theta, phi)
T p(T theta, T phi) {
T fx = 0;
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;
}
|
487a9b49
David Mayerich
added the ability...
|
268
|
}
|
01707489
David Mayerich
implemented a new...
|
269
|
return fx;
|
487a9b49
David Mayerich
added the ability...
|
270
271
|
}
|
01707489
David Mayerich
implemented a new...
|
272
273
274
275
276
277
|
/// Returns the derivative of the spherical function with respect to theta
/// return value is in cartesian coordinates
vec3<T> 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);
|
487a9b49
David Mayerich
added the ability...
|
278
|
|
01707489
David Mayerich
implemented a new...
|
279
280
281
|
vec3<T> s(r, theta, phi); //get the spherical coordinate position for all three points
vec3<T> st(rt, theta + d, phi);
//vec3<double> sp(rp, theta, phi + d);
|
487a9b49
David Mayerich
added the ability...
|
282
|
|
01707489
David Mayerich
implemented a new...
|
283
284
285
|
vec3<T> c = s.sph2cart();
vec3<T> ct = st.sph2cart();
//vec3<double> cp = sp.sph2cart();
|
487a9b49
David Mayerich
added the ability...
|
286
|
|
01707489
David Mayerich
implemented a new...
|
287
288
289
|
vec3<T> dt = (ct - c)/d; //calculate the derivative
return dt;
}
|
487a9b49
David Mayerich
added the ability...
|
290
|
|
01707489
David Mayerich
implemented a new...
|
291
292
293
294
295
296
|
/// Returns the derivative of the spherical function with respect to phi
/// return value is in cartesian coordinates
vec3<T> 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);
|
487a9b49
David Mayerich
added the ability...
|
297
|
|
01707489
David Mayerich
implemented a new...
|
298
299
300
|
vec3<T> s(r, theta, phi); //get the spherical coordinate position for all three points
//vec3<double> st(rt, theta + d, phi);
vec3<T> sp(rp, theta, phi + d);
|
487a9b49
David Mayerich
added the ability...
|
301
|
|
01707489
David Mayerich
implemented a new...
|
302
303
304
|
vec3<T> c = s.sph2cart();
//vec3<double> ct = st.sph2cart();
vec3<T> cp = sp.sph2cart();
|
487a9b49
David Mayerich
added the ability...
|
305
|
|
01707489
David Mayerich
implemented a new...
|
306
307
|
vec3<T> dp = (cp - c) / d; //calculate the derivative
return dp;
|
487a9b49
David Mayerich
added the ability...
|
308
|
}
|
01707489
David Mayerich
implemented a new...
|
309
310
311
312
313
314
|
/// 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);
|
424c1475
David Mayerich
added a grid samp...
|
315
|
}
|
01707489
David Mayerich
implemented a new...
|
316
317
318
319
320
321
322
323
324
|
/// 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);
}
}
|
761ebaa9
David Mayerich
spherical harmoni...
|
325
|
}
|
761ebaa9
David Mayerich
spherical harmoni...
|
326
|
|
01707489
David Mayerich
implemented a new...
|
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
|
/// 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
|
487a9b49
David Mayerich
added the ability...
|
347
348
349
350
351
352
353
|
}
|
01707489
David Mayerich
implemented a new...
|
354
|
#endif
|