sphericalHarmonics.hpp 11.8 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

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)];
19
template <typename Tscalar,typename Tcoef> Tscalar SphericalHarmonics<Tscalar,Tcoef>::F_tab[32][(max_resolution+1)*(max_resolution+1)];
20 21 22 23 24 25 26 27 28


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

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

/*************************************************************************
evaluation
**************************************************************************/
78

79
template <typename Tscalar,typename Tcoef>
80
void SphericalHarmonics<Tscalar,Tcoef>::set_eval_direction (Tscalar theta, Tscalar phi, unsigned int threadId)
81
{
82 83 84
	assert(threadId < 32);
	compute_P_tab(std::cos(theta), threadId);
	compute_y_tab(phi, threadId);
85 86 87
}

template <typename Tscalar,typename Tcoef>
88
void SphericalHarmonics<Tscalar,Tcoef>::set_eval_direction (Tscalar x, Tscalar y, Tscalar z, unsigned int threadId)
89
{
90 91
	assert(threadId < 32);
	compute_P_tab(z, threadId);
92 93

	Tscalar phi (0);
94 95
	if ((x*x + y*y) > 0.0)
		phi = atan2(y, x);
96

97
	compute_y_tab(phi, threadId);
98 99 100
}

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

template <typename Tscalar,typename Tcoef>
113
Tcoef SphericalHarmonics<Tscalar,Tcoef>::evaluate_at (Tscalar theta, Tscalar phi, unsigned int threadId) const
114
{
115 116 117
	assert(threadId < 32);
	set_eval_direction(theta, phi, threadId);
	return evaluate(threadId);
118 119 120
}

template <typename Tscalar,typename Tcoef>
121
Tcoef SphericalHarmonics<Tscalar,Tcoef>::evaluate_at (Tscalar x, Tscalar y, Tscalar z, unsigned int threadId) const
122
{
123 124 125
	assert(threadId < 32);
	set_eval_direction(x, y, z, threadId);
	return evaluate(threadId);
126 127 128 129 130
}

template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::init_K_tab ()
{
131
	for(int l = 0; l <= resolution; l++)
132 133
	{
		// recursive computation of the squares
134
		K_tab[index(l,0)] = Tscalar((2*l+1) / (4*M_PI));
135
		for (int m = 1; m <= l; m++)
136 137 138
			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)]);
139
		for (int m = 1; m <= l; m++)
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>
158
void SphericalHarmonics<Tscalar,Tcoef>::compute_P_tab (Tscalar t, unsigned int threadId)
159 160 161
{
//	if (t<0) {t=-t*t;} else {t=t*t;} // for plotting only : expand the param near equator

162
	F_tab[threadId][index(0,0)] = 1;
163
	for (int l = 1; l <= resolution; l++)
164
	{
165 166
		F_tab[threadId][index(l,l)] = (1-2*l) * sqrt(1-t*t) * F_tab[threadId][index(l-1,l-1)];  // first diago
		F_tab[threadId][index(l,l-1)] = t * (2*l-1) * F_tab[threadId][index(l-1,l-1)];// second diago
167
		for (int m = 0; m <= l-2; m++)
168
		{// remaining of the line under the 2 diago
169
			F_tab[threadId][index(l,m)] = t * (2*l-1) / (float) (l-m) * F_tab[threadId][index(l-1,m)] - (l+m-1) / (float) (l-m) * F_tab[threadId][index(l-2,m)];
170 171 172 173 174
		}
	}
}

template <typename Tscalar,typename Tcoef>
175
void SphericalHarmonics<Tscalar,Tcoef>::compute_y_tab (Tscalar phi, unsigned int threadId)
176
{
177
	for (int l = 0; l <= resolution; l++)
178
	{
179
		F_tab[threadId][index(l,0)] *= K_tab[index(l,0)]; // remove for plotting
180 181
	}

182
	for (int m = 1; m <= resolution; m++)
183
	{
184 185
		Tscalar cos_m_phi = std::cos ( m * phi );
		Tscalar sin_m_phi = std::sin ( m * phi );
186

187
		for (int l = m; l <= resolution; l++)
188
		{
189 190 191 192
			F_tab[threadId][index(l,m)] *= std::sqrt(2.0); // remove for plotting
			F_tab[threadId][index(l,m)] *= K_tab[index(l,m)]; // remove for plotting
			F_tab[threadId][index(l,-m)] = F_tab[threadId][index(l,m)] * sin_m_phi ; // store the values for -m<0 in the upper triangle
			F_tab[threadId][index(l,m)] *= cos_m_phi;
193 194 195 196 197 198 199 200
		}
	}

}

/*************************************************************************
I/O
**************************************************************************/
201

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

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

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

template <typename Tscalar,typename Tcoef>
226
void SphericalHarmonics<Tscalar,Tcoef>::operator += (const SphericalHarmonics<Tscalar,Tcoef>& sh)
227
{
228
	for (int i = 0; i < nb_coefs; i++)
229 230 231 232
		this->coefs[i] += sh.coefs[i];
}

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

247 248 249 250 251 252 253 254
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;
}

255 256 257
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::operator *= (Tscalar s)
{
258
	for (int i = 0; i < nb_coefs; i++)
259 260 261
		this->coefs[i] *= s;
}

262 263 264 265 266 267 268 269
template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef> SphericalHarmonics<Tscalar,Tcoef>::operator * (Tscalar s) const
{
	SphericalHarmonics<Tscalar,Tcoef> res(*this);
	res *= s;
	return res;
}

270 271 272
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::operator /= (Tscalar s)
{
273
	for (int i = 0; i < nb_coefs; i++)
274 275 276
		this->coefs[i] /= s;
}

277 278 279 280 281 282 283 284
template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef> SphericalHarmonics<Tscalar,Tcoef>::operator / (Tscalar s) const
{
	SphericalHarmonics<Tscalar,Tcoef> res(*this);
	res /= s;
	return res;
}

285 286 287 288 289 290
/*************************************************************************
fitting
**************************************************************************/

template <typename Tscalar,typename Tcoef>
template <typename Tdirection, typename Tchannel>
291 292 293 294
void SphericalHarmonics<Tscalar,Tcoef>::fit_to_data(
	int n,
	Tdirection* t_theta, Tdirection* t_phi,
	Tchannel* t_R, Tchannel* t_G, Tchannel* t_B,
295 296
	double lambda,
	unsigned int threadId)
297
{
298
	Eigen::MatrixXd mM (nb_coefs,n); // matrix with basis function values, evaluated for all directions
299
	// compute mM
300
	for (int p = 0; p < n; ++p)
301
	{
302
		set_eval_direction(t_theta[p], t_phi[p], threadId);
303
		for (int i = 0; i < nb_coefs; ++i)
304
		{
305
			mM(i,p) = F_tab[threadId][i];
306 307
		}
	}
308
	fit_to_data(n, mM, t_R, t_G, t_B, lambda);
309 310 311 312
}

template <typename Tscalar,typename Tcoef>
template <typename Tdirection, typename Tchannel>
313 314 315 316
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,
317 318
	double lambda,
	unsigned int threadId)
319
{
320
	Eigen::MatrixXd mM (nb_coefs,n); // matrix with basis function values, evaluated for all directions
321 322 323
	// compute mM
	for (int p=0; p<n; ++p)
	{
324
		set_eval_direction(t_x[p],t_y[p],t_z[p], threadId);
325 326
		for (int i=0; i<nb_coefs; ++i)
		{
327
			mM(i,p) = F_tab[threadId][i];
328 329
		}
	}
330
	fit_to_data(n, mM, t_R, t_G, t_B, lambda);
331 332 333 334
}

template <typename Tscalar,typename Tcoef>
template <typename Tchannel>
335 336 337 338 339
void SphericalHarmonics<Tscalar,Tcoef>::fit_to_data(
	int n,
	Eigen::MatrixXd& mM,
	Tchannel* t_R, Tchannel* t_G, Tchannel* t_B,
	double lambda)
340 341 342 343 344 345
{
	// 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
346 347 348
	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
349 350

	// compute mA
351
	for (int i = 0; i < nb_coefs; ++i)
352
	{
353
		for (int j = 0; j < nb_coefs; ++j)
354 355
		{
			mA(i,j) = 0;
356
			for (int p = 0; p < n; ++p)
357 358 359 360 361 362 363
			{
				mA(i,j) += mM(i,p) * mM(j,p);
			}
			mA(i,j) *= (1.0-lambda) / n;
		}
	}

364
	for (int l = 0; l <= resolution; ++l)
365
	{
366
		for (int m =- l; m <= l; ++m)
367 368 369 370 371 372 373
		{
			int i = index(l,m);
			mA(i,i) += lambda * l * (l+1) / (4.0*M_PI);
		}
	}

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

411
} // namespace Utils
412

413
} // namespace CGoGN