sphericalHarmonics.hpp 11.1 KB
Newer Older
1
2
3
4
5
6
7
/*
 * sphericalHarmonics.hpp
 *
 *  Created on: Oct 2, 2013
 *      Author: sauvage
 */

8
9
10
11
12
namespace CGoGN
{

namespace Utils
{
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

template <typename Tscalar,typename Tcoef> int SphericalHarmonics<Tscalar,Tcoef>::resolution = -1;
template <typename Tscalar,typename Tcoef> int SphericalHarmonics<Tscalar,Tcoef>::nb_coefs = -1;
template <typename Tscalar,typename Tcoef> unsigned long SphericalHarmonics<Tscalar,Tcoef>::cpt_instances = 0;

template <typename Tscalar,typename Tcoef> Tscalar SphericalHarmonics<Tscalar,Tcoef>::K_tab[(max_resolution+1)*(max_resolution+1)];
template <typename Tscalar,typename Tcoef> Tscalar SphericalHarmonics<Tscalar,Tcoef>::F_tab[(max_resolution+1)*(max_resolution+1)];


/*************************************************************************
construction, destruction and initialization
**************************************************************************/

template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef>::SphericalHarmonics()
{
29
30
	assert ( (nb_coefs > 0) || !" maybe you forgot to call set_level()");
	coefs = new Tcoef[nb_coefs];
31
	++ cpt_instances;
32
33
	for (int i = 0; i < nb_coefs; i++)
		coefs[i] = Tcoef (0);
34
35
36
37
38
}

template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef>::SphericalHarmonics(SphericalHarmonics const & r)
{
39
40
41
42
	assert ( (nb_coefs > 0) || !" maybe you forgot to call set_level()");
	coefs = new Tcoef[nb_coefs];
	++cpt_instances;
	for (int i = 0; i < nb_coefs; i++)
43
44
45
46
47
48
		coefs[i] = r.coefs[i];
}

template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef>::~SphericalHarmonics()
{
49
50
	delete[] coefs;
	--cpt_instances;
51
52
53
54
55
}

template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::set_level(int res_level)
{
56
57
	assert(res_level >= 0 && res_level < max_resolution);
	assert(cpt_instances == 0);
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
	resolution = res_level;
	nb_coefs = (resolution + 1) * (resolution + 1);
	init_K_tab();
}

template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::set_nb_coefs(int nbc)
{
	assert(nbc > 0);
	int sq = ceil(sqrt(nbc)) ;
	assert(sq*sq == nbc || !"Number of coefs does not fill the last level") ;
	set_level(sq-1) ;
}

/*************************************************************************
evaluation
**************************************************************************/
75

76
77
78
79
80
81
82
83
84
85
86
87
88
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::set_eval_direction (Tscalar theta, Tscalar phi)
{
	compute_P_tab(cos(theta));
	compute_y_tab(phi);
}

template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::set_eval_direction (Tscalar x, Tscalar y, Tscalar z)
{
	compute_P_tab(z);

	Tscalar phi (0);
89
90
	if ((x*x + y*y) > 0.0)
		phi = atan2(y, x);
91
92
93
94
95
96
97
98

	compute_y_tab(phi);
}

template <typename Tscalar,typename Tcoef>
Tcoef SphericalHarmonics<Tscalar,Tcoef>::evaluate () const
{
	Tcoef r (0); // (0.0,0.0,0.0); //  TODO : use Tcoef (0)
99
	for (int i = 0; i < nb_coefs; i++)
100
101
102
103
104
105
106
107
108
	{
		r += coefs[i] * F_tab[i];
	}
	return r;
}

template <typename Tscalar,typename Tcoef>
Tcoef SphericalHarmonics<Tscalar,Tcoef>::evaluate_at (Tscalar theta, Tscalar phi) const
{
109
	set_eval_direction(theta, phi);
110
111
112
113
114
115
	return evaluate();
}

template <typename Tscalar,typename Tcoef>
Tcoef SphericalHarmonics<Tscalar,Tcoef>::evaluate_at (Tscalar x, Tscalar y, Tscalar z) const
{
116
	set_eval_direction(x, y, z);
117
118
119
120
121
122
	return evaluate();
}

template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::init_K_tab ()
{
123
	for(int l = 0; l <= resolution; l++)
124
125
126
	{
		// recursive computation of the squares
		K_tab[index(l,0)] = (2*l+1) / (4*M_PI);
127
		for (int m = 1; m <= l; m++)
128
129
130
			K_tab[index(l,m)] = K_tab[index(l,m-1)] / (l-m+1) / (l+m);
		// square root + symmetry
		K_tab[index(l,0)] = sqrt(K_tab[index(l,0)]);
131
		for (int m = 1; m <= l; m++)
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
		{
			K_tab[index(l,m)] = sqrt(K_tab[index(l,m)]);
			K_tab[index(l,-m)] = K_tab[index(l,m)];
		}
	}
}

/* obsolete : was used for shaders
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::copy_K_tab (Tscalar tab[])
{
	assert ( (nb_coefs>0) || !" maybe you forgot to call set_level()");
	for (unsigned int i = 0 ; i < nb_coefs ; ++i)
		tab[i] = K_tab[i] ;
}
*/

template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::compute_P_tab (Tscalar t)
{
//	if (t<0) {t=-t*t;} else {t=t*t;} // for plotting only : expand the param near equator

	F_tab[index(0,0)] = 1;
155
	for (int l = 1; l <= resolution; l++)
156
157
158
	{
		F_tab[index(l,l)] = (1-2*l) * sqrt(1-t*t) * F_tab[index(l-1,l-1)];  // first diago
		F_tab[index(l,l-1)] = t * (2*l-1) * F_tab[index(l-1,l-1)];// second diago
159
		for (int m = 0; m <= l-2; m++)
160
161
162
163
164
165
166
167
168
		{// remaining of the line under the 2 diago
			F_tab[index(l,m)] = t * (2*l-1) / (float) (l-m) * F_tab[index(l-1,m)] - (l+m-1) / (float) (l-m) * F_tab[index(l-2,m)];
		}
	}
}

template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::compute_y_tab (Tscalar phi)
{
169
	for (int l = 0; l <= resolution; l++)
170
171
172
173
	{
		F_tab[index(l,0)] *= K_tab[index(l,0)]; // remove for plotting
	}

174
	for (int m = 1; m <= resolution; m++)
175
176
177
178
	{
		Tscalar cos_m_phi = cos ( m * phi );
		Tscalar sin_m_phi = sin ( m * phi );

179
		for (int l = m; l <= resolution; l++)
180
181
182
183
184
185
186
187
188
189
190
191
192
		{
			F_tab[index(l,m)] *= sqrt(2.0); // remove for plotting
			F_tab[index(l,m)] *= K_tab[index(l,m)]; // remove for plotting
			F_tab[index(l,-m)] = F_tab[index(l,m)] * sin_m_phi ; // store the values for -m<0 in the upper triangle
			F_tab[index(l,m)] *= cos_m_phi;
		}
	}

}

/*************************************************************************
I/O
**************************************************************************/
193

194
195
196
template <typename Tscalar,typename Tcoef>
std::ostream & operator << (std::ostream & os, const SphericalHarmonics<Tscalar,Tcoef> & sh)
{
197
	for (int l = 0; l <= sh.resolution; l++)
198
	{
199
		for (int m = -l; m <= l; m++)
200
201
202
203
204
205
206
207
208
209
210
			os << sh.get_coef(l,m) << "\t";
		os << std::endl;
	}
	return os;
}

/*************************************************************************
operators
**************************************************************************/

template <typename Tscalar,typename Tcoef>
211
void SphericalHarmonics<Tscalar,Tcoef>::operator = (const SphericalHarmonics<Tscalar,Tcoef>& sh)
212
{
213
	for (int i = 0; i < nb_coefs; i++)
214
215
216
217
		this->coefs[i] = sh.coefs[i];
}

template <typename Tscalar,typename Tcoef>
218
void SphericalHarmonics<Tscalar,Tcoef>::operator += (const SphericalHarmonics<Tscalar,Tcoef>& sh)
219
{
220
	for (int i = 0; i < nb_coefs; i++)
221
222
223
224
		this->coefs[i] += sh.coefs[i];
}

template <typename Tscalar,typename Tcoef>
225
SphericalHarmonics<Tscalar,Tcoef> SphericalHarmonics<Tscalar,Tcoef>::operator + (const SphericalHarmonics<Tscalar,Tcoef>& sh) const
226
{
227
228
229
230
231
232
233
234
235
	SphericalHarmonics<Tscalar,Tcoef> res(*this);
	res += sh;
	return res;
}

template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::operator -= (const SphericalHarmonics<Tscalar,Tcoef>& sh)
{
	for (int i = 0; i < nb_coefs; i++)
236
237
238
		this->coefs[i] -= sh.coefs[i];
}

239
240
241
242
243
244
245
246
template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef> SphericalHarmonics<Tscalar,Tcoef>::operator - (const SphericalHarmonics<Tscalar,Tcoef>& sh) const
{
	SphericalHarmonics<Tscalar,Tcoef> res(*this);
	res -= sh;
	return res;
}

247
248
249
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::operator *= (Tscalar s)
{
250
	for (int i = 0; i < nb_coefs; i++)
251
252
253
		this->coefs[i] *= s;
}

254
255
256
257
258
259
260
261
template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef> SphericalHarmonics<Tscalar,Tcoef>::operator * (Tscalar s) const
{
	SphericalHarmonics<Tscalar,Tcoef> res(*this);
	res *= s;
	return res;
}

262
263
264
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::operator /= (Tscalar s)
{
265
	for (int i = 0; i < nb_coefs; i++)
266
267
268
		this->coefs[i] /= s;
}

269
270
271
272
273
274
275
276
template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef> SphericalHarmonics<Tscalar,Tcoef>::operator / (Tscalar s) const
{
	SphericalHarmonics<Tscalar,Tcoef> res(*this);
	res /= s;
	return res;
}

277
278
279
280
281
282
/*************************************************************************
fitting
**************************************************************************/

template <typename Tscalar,typename Tcoef>
template <typename Tdirection, typename Tchannel>
283
284
285
286
287
void SphericalHarmonics<Tscalar,Tcoef>::fit_to_data(
	int n,
	Tdirection* t_theta, Tdirection* t_phi,
	Tchannel* t_R, Tchannel* t_G, Tchannel* t_B,
	double lambda)
288
{
289
	Eigen::MatrixXd mM (nb_coefs,n); // matrix with basis function values, evaluated for all directions
290
	// compute mM
291
	for (int p = 0; p < n; ++p)
292
	{
293
294
		set_eval_direction(t_theta[p], t_phi[p]);
		for (int i = 0; i < nb_coefs; ++i)
295
296
297
298
		{
			mM(i,p) = F_tab[i];
		}
	}
299
	fit_to_data(n, mM, t_R, t_G, t_B, lambda);
300
301
302
303
}

template <typename Tscalar,typename Tcoef>
template <typename Tdirection, typename Tchannel>
304
305
306
307
308
void SphericalHarmonics<Tscalar,Tcoef>::fit_to_data(
	int n,
	Tdirection* t_x, Tdirection* t_y, Tdirection* t_z,
	Tchannel* t_R, Tchannel* t_G, Tchannel* t_B,
	double lambda)
309
{
310
	Eigen::MatrixXd mM (nb_coefs,n); // matrix with basis function values, evaluated for all directions
311
312
313
314
315
316
317
318
319
	// compute mM
	for (int p=0; p<n; ++p)
	{
		set_eval_direction(t_x[p],t_y[p],t_z[p]);
		for (int i=0; i<nb_coefs; ++i)
		{
			mM(i,p) = F_tab[i];
		}
	}
320
	fit_to_data(n, mM, t_R, t_G, t_B, lambda);
321
322
323
324
}

template <typename Tscalar,typename Tcoef>
template <typename Tchannel>
325
326
327
328
329
void SphericalHarmonics<Tscalar,Tcoef>::fit_to_data(
	int n,
	Eigen::MatrixXd& mM,
	Tchannel* t_R, Tchannel* t_G, Tchannel* t_B,
	double lambda)
330
331
332
333
334
335
{
	// fits the data t_R, t_G and t_B, according to our 2013 CGF paper
	// mM contains basis function values, (already) evaluated for all input directions
	// works only for 3 channels

	// allocate the memory
336
337
338
	Eigen::MatrixXd mA (nb_coefs, nb_coefs); // matrix A in linear system AC=B
	Eigen::MatrixXd mB (nb_coefs, 3); // matrix B in linear system AC=B : contains [t_R, t_G, t_B]
	Eigen::MatrixXd mC (nb_coefs, 3); // matrix C (solution) in linear system AC=B : contains the RGB coefs of the resulting SH
339
340

	// compute mA
341
	for (int i = 0; i < nb_coefs; ++i)
342
	{
343
		for (int j = 0; j < nb_coefs; ++j)
344
345
		{
			mA(i,j) = 0;
346
			for (int p = 0; p < n; ++p)
347
348
349
350
351
352
353
			{
				mA(i,j) += mM(i,p) * mM(j,p);
			}
			mA(i,j) *= (1.0-lambda) / n;
		}
	}

354
	for (int l = 0; l <= resolution; ++l)
355
	{
356
		for (int m =- l; m <= l; ++m)
357
358
359
360
361
362
363
		{
			int i = index(l,m);
			mA(i,i) += lambda * l * (l+1) / (4.0*M_PI);
		}
	}

	// compute mB
364
	for (int i = 0; i < nb_coefs; ++i)
365
366
367
368
	{
		mB(i,0) = 0.0;
		mB(i,1) = 0.0;
		mB(i,2) = 0.0;
369
		for (int p = 0; p < n; ++p)
370
		{
371
372
373
			mB(i,0) += mM(i,p) * t_R[p];
			mB(i,1) += mM(i,p) * t_G[p];
			mB(i,2) += mM(i,p) * t_B[p];
374
375
376
377
378
379
380
		}
		mB(i,0) *= (1.0-lambda) / n;
		mB(i,1) *= (1.0-lambda) / n;
		mB(i,2) *= (1.0-lambda) / n;
	}

	// solve the system with LDLT decomposition
381
	Eigen::LDLT<Eigen::MatrixXd> solver (mA);
382
383
384
385
386
387
388
389
390
391
392
	mC = solver.solve(mB);

//	std::cout << "phi = " << mM << std::endl;
//	std::cout << mA << std::endl;
//	std::cout << "*" << std::endl;
//	std::cout << mC << std::endl;
//	std::cout << "=" << std::endl;
//	std::cout << mB << std::endl;

	// store result in the SH
	// it is assumed that Tcoef is VEC3 actually
393
	for (int i = 0; i < nb_coefs; ++i)
394
395
396
397
398
399
400
	{
		get_coef(i)[0] = mC(i,0);
		get_coef(i)[1] = mC(i,1);
		get_coef(i)[2] = mC(i,2);
	}
}

401
} // namespace Utils
402

403
} // namespace CGoGN