Création d'un compte pour un collaborateur extérieur au laboratoire depuis l'intranet ICube : https://intranet.icube.unistra.fr/fr/labs/member/profile

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
58
59
60
61
62
63
	if (res_level != resolution)
	{
		assert(res_level >= 0 && res_level < max_resolution);
		assert(cpt_instances == 0);
		resolution = res_level;
		nb_coefs = (resolution + 1) * (resolution + 1);
		init_K_tab();
	}
64
65
66
67
68
69
70
71
72
73
74
75
76
77
}

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
**************************************************************************/
78

79
80
81
82
83
84
85
86
87
88
89
90
91
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);
92
93
	if ((x*x + y*y) > 0.0)
		phi = atan2(y, x);
94
95
96
97
98
99
100
101

	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)
102
	for (int i = 0; i < nb_coefs; i++)
103
104
105
106
107
108
109
110
111
	{
		r += coefs[i] * F_tab[i];
	}
	return r;
}

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

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

template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::init_K_tab ()
{
126
	for(int l = 0; l <= resolution; l++)
127
128
129
	{
		// recursive computation of the squares
		K_tab[index(l,0)] = (2*l+1) / (4*M_PI);
130
		for (int m = 1; m <= l; m++)
131
132
133
			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)]);
134
		for (int m = 1; m <= l; m++)
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
		{
			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;
158
	for (int l = 1; l <= resolution; l++)
159
160
161
	{
		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
162
		for (int m = 0; m <= l-2; m++)
163
164
165
166
167
168
169
170
171
		{// 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)
{
172
	for (int l = 0; l <= resolution; l++)
173
174
175
176
	{
		F_tab[index(l,0)] *= K_tab[index(l,0)]; // remove for plotting
	}

177
	for (int m = 1; m <= resolution; m++)
178
179
180
181
	{
		Tscalar cos_m_phi = cos ( m * phi );
		Tscalar sin_m_phi = sin ( m * phi );

182
		for (int l = m; l <= resolution; l++)
183
184
185
186
187
188
189
190
191
192
193
194
195
		{
			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
**************************************************************************/
196

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

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

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

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

template <typename Tscalar,typename Tcoef>
228
SphericalHarmonics<Tscalar,Tcoef> SphericalHarmonics<Tscalar,Tcoef>::operator + (const SphericalHarmonics<Tscalar,Tcoef>& sh) const
229
{
230
231
232
233
234
235
236
237
238
	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++)
239
240
241
		this->coefs[i] -= sh.coefs[i];
}

242
243
244
245
246
247
248
249
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;
}

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

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

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

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

280
281
282
283
284
285
/*************************************************************************
fitting
**************************************************************************/

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

template <typename Tscalar,typename Tcoef>
template <typename Tdirection, typename Tchannel>
307
308
309
310
311
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)
312
{
313
	Eigen::MatrixXd mM (nb_coefs,n); // matrix with basis function values, evaluated for all directions
314
315
316
317
318
319
320
321
322
	// 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];
		}
	}
323
	fit_to_data(n, mM, t_R, t_G, t_B, lambda);
324
325
326
327
}

template <typename Tscalar,typename Tcoef>
template <typename Tchannel>
328
329
330
331
332
void SphericalHarmonics<Tscalar,Tcoef>::fit_to_data(
	int n,
	Eigen::MatrixXd& mM,
	Tchannel* t_R, Tchannel* t_G, Tchannel* t_B,
	double lambda)
333
334
335
336
337
338
{
	// 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
339
340
341
	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
342
343

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

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

	// compute mB
367
	for (int i = 0; i < nb_coefs; ++i)
368
369
370
371
	{
		mB(i,0) = 0.0;
		mB(i,1) = 0.0;
		mB(i,2) = 0.0;
372
		for (int p = 0; p < n; ++p)
373
		{
374
375
376
			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];
377
378
379
380
381
382
383
		}
		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
384
	Eigen::LDLT<Eigen::MatrixXd> solver (mA);
385
386
387
388
389
390
391
392
393
394
395
	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
396
	for (int i = 0; i < nb_coefs; ++i)
397
398
399
400
401
402
403
	{
		get_coef(i)[0] = mC(i,0);
		get_coef(i)[1] = mC(i,1);
		get_coef(i)[2] = mC(i,2);
	}
}

404
} // namespace Utils
405

406
} // namespace CGoGN