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