/*

Copyright (c) 2004 Jan Kautz

This software is provided 'as-is', without any express or implied
warranty. In no event will the authors be held liable for any
damages arising from the use of this software.

Permission is granted to anyone to use this software for any
*non-commerical* purpose. 

*/

// WARNING: this is only a code-snippet to show the creation and use of SH gradients, the code does not compile!

// ------------------------------------------------
// Code to create gradients
// ------------------------------------------------

{
  {
	// clear gradients (basically an array of coefficients (RGB each))
	gradientshX.clear();
	gradientshY.clear();
	gradientshZ.clear();

	// room for gradients
	float dYdx[9], dYdy[9], dYdz[9];

	// setup SH
	static SimpleSphericalHarmonics s(Lmax);
	static float ylm[25]; // will store SH evaluated for a directions

	// iterate over environment maps
	CubeEnvironmentMapIterator<unsigned char,3> it(readCubeMap);
	CubeEnvironmentMapIterator<float,1> dit(depthCubeMap);
	//CubeEnvironmentMapIterator<float,3> nit(normalCubeMap); // in case normal is used
	for( it.begin(), dit.begin(); !it.end(); ++it, ++dit ) {
	  // get direction of current env-map texel 
	  Vec3f dir = it.getDirection();
	  float r = dit.getValue()[0]; // get depth-value for that texel
	  Vec3f rdir = dir * r; // 
	  //Vec3f n_s = nit.getValue(); // get normal (accurate version, but requires cube map with normals in it)
	  Vec3f n_s = dir; // simplified version: assumes all surfaces face towards center of envmap
	 
	  float nx = n_s[0], ny = n_s[1], nz = n_s[2];
	  float x = rdir[0], y = rdir[1], z = rdir[2];
	  float ds = it.getDiffSolidAngle(); // get solid angle for current cubemap texel
	  Vec3f I_s = it.getValue()/255.f;   // get intensity (RGB)

	  // compute ylmArray
	  s.computeYlmArray( dir, ylm ); // fill ylm[] with SH coefficients evaluated for direction dir
	  computedYdx( rdir, dYdx );     // compute gradients in x,y,z
	  computedYdy( rdir, dYdy );
	  computedYdz( rdir, dYdz );

	  const int R = 0, G = 1, B = 2;

	  if( gradientMode == 0 ) {
		// simplified version, no normal -> computation is easier and faster

		for( int i = 0; i < 9; ++i ) {
		  //grad(lambda_i) += y_i(s/||s||)*I_s*ds*r^2*grad(1/||s||^2) + I_s*ds * grad(y_i(s/||s||))
		  (*gradientshX.getCoefficients(i))[R] += ylm[i]*I_s[R]*ds*r*r* (-2*x)/(pow(x*x+y*y+z*z,2.0f)) +
			I_s[R]*ds * dYdx[i];
		  (*gradientshX.getCoefficients(i))[G] += ylm[i]*I_s[G]*ds*r*r* (-2*x)/(pow(x*x+y*y+z*z,2.0f)) +
			I_s[G]*ds * dYdx[i];
		  (*gradientshX.getCoefficients(i))[B] += ylm[i]*I_s[B]*ds*r*r* (-2*x)/(pow(x*x+y*y+z*z,2.0f)) +
			I_s[B]*ds * dYdx[i];

		  (*gradientshY.getCoefficients(i))[R] += ylm[i]*I_s[R]*ds*r*r* (-2*y)/(pow(x*x+y*y+z*z,2.0f)) +
			I_s[R]*ds * dYdy[i];
		  (*gradientshY.getCoefficients(i))[G] += ylm[i]*I_s[G]*ds*r*r* (-2*y)/(pow(x*x+y*y+z*z,2.0f)) +
			I_s[G]*ds * dYdy[i];
		  (*gradientshY.getCoefficients(i))[B] += ylm[i]*I_s[B]*ds*r*r* (-2*y)/(pow(x*x+y*y+z*z,2.0f)) +
			I_s[B]*ds * dYdy[i];

		  (*gradientshZ.getCoefficients(i))[R] += ylm[i]*I_s[R]*ds*r*r* (-2*z)/(pow(x*x+y*y+z*z,2.0f)) +
			I_s[R]*ds * dYdz[i];
		  (*gradientshZ.getCoefficients(i))[G] += ylm[i]*I_s[G]*ds*r*r* (-2*z)/(pow(x*x+y*y+z*z,2.0f)) +
			I_s[G]*ds * dYdz[i];
		  (*gradientshZ.getCoefficients(i))[B] += ylm[i]*I_s[B]*ds*r*r* (-2*z)/(pow(x*x+y*y+z*z,2.0f)) +
			I_s[B]*ds * dYdz[i];
		}
	  } else {
		// use normals

		for( int i = 0; i < 9; ++i ) {
		  //grad(lambda_i) += y_i(s/||s||)*I_s*ds*r^2/(n_s dot sN) * grad((n dot s/||s||)/||s||^2) + I_s*ds * grad(y_i(s/||s||))
		  (*gradientshX.getCoefficients(i))[R] += ylm[i]*I_s[R]*ds*r*r/(n_s|dir) * 
			(nx/sqrt(pow(x*x+y*y+z*z,3.0f)) - (3*x)*(nx*x+ny*y+nz*z)/sqrt(pow(x*x+y*y+z*z,5.0f))) +
			I_s[R]*ds * dYdx[i];
		  (*gradientshX.getCoefficients(i))[G] += ylm[i]*I_s[G]*ds*r*r/(n_s|dir) * 
			(nx/sqrt(pow(x*x+y*y+z*z,3.0f)) - (3*x)*(nx*x+ny*y+nz*z)/sqrt(pow(x*x+y*y+z*z,5.0f))) +
			I_s[G]*ds * dYdx[i];
		  (*gradientshX.getCoefficients(i))[B] += ylm[i]*I_s[B]*ds*r*r/(n_s|dir) * 
			(nx/sqrt(pow(x*x+y*y+z*z,3.0f)) - (3*x)*(nx*x+ny*y+nz*z)/sqrt(pow(x*x+y*y+z*z,5.0f))) +
			I_s[B]*ds * dYdx[i];

		  (*gradientshY.getCoefficients(i))[R] += ylm[i]*I_s[R]*ds*r*r/(n_s|dir) * 
			(ny/sqrt(pow(x*x+y*y+z*z,3.0f)) - (3*y)*(nx*x+ny*y+nz*z)/sqrt(pow(x*x+y*y+z*z,5.0f))) +
			I_s[R]*ds * dYdy[i];
		  (*gradientshY.getCoefficients(i))[G] += ylm[i]*I_s[G]*ds*r*r/(n_s|dir) * 
			(ny/sqrt(pow(x*x+y*y+z*z,3.0f)) - (3*y)*(nx*x+ny*y+nz*z)/sqrt(pow(x*x+y*y+z*z,5.0f))) +
			I_s[G]*ds * dYdy[i];
		  (*gradientshY.getCoefficients(i))[B] += ylm[i]*I_s[B]*ds*r*r/(n_s|dir) * 
			(ny/sqrt(pow(x*x+y*y+z*z,3.0f)) - (3*y)*(nx*x+ny*y+nz*z)/sqrt(pow(x*x+y*y+z*z,5.0f))) +
			I_s[B]*ds * dYdy[i];

		  (*gradientshZ.getCoefficients(i))[R] += ylm[i]*I_s[R]*ds*r*r/(n_s|dir) * 
			(nz/sqrt(pow(x*x+y*y+z*z,3.0f)) - (3*z)*(nx*x+ny*y+nz*z)/sqrt(pow(x*x+y*y+z*z,5.0f))) +
			I_s[R]*ds * dYdz[i];
		  (*gradientshZ.getCoefficients(i))[G] += ylm[i]*I_s[G]*ds*r*r/(n_s|dir) * 
			(nz/sqrt(pow(x*x+y*y+z*z,3.0f)) - (3*z)*(nx*x+ny*y+nz*z)/sqrt(pow(x*x+y*y+z*z,5.0f))) +
			I_s[G]*ds * dYdz[i];
		  (*gradientshZ.getCoefficients(i))[B] += ylm[i]*I_s[B]*ds*r*r/(n_s|dir) * 
			(nz/sqrt(pow(x*x+y*y+z*z,3.0f)) - (3*z)*(nx*x+ny*y+nz*z)/sqrt(pow(x*x+y*y+z*z,5.0f))) +
			I_s[B]*ds * dYdz[i];
		}
	  }
	}

	// now invert the gradients (computed them in the wrong direction)
	for( i = 0; i < 9; ++i ) {
	  gradientshX.getCoefficients()[i] = -gradientshX.getCoefficients()[i];
	  gradientshY.getCoefficients()[i] = -gradientshY.getCoefficients()[i];
	  gradientshZ.getCoefficients()[i] = -gradientshZ.getCoefficients()[i];	 
	}
  }
}

// ------------------------------------------------
// Actual gradient computation
// ------------------------------------------------

void
GradientSH::computedYdx( Vec3f dir, float *dYdx ) {
  float t1, t2, t3, t4, t5, t7, t19, t23, t24;
  float x = dir[0];
  float y = dir[1];
  float z = dir[2];

  t1 = x*x;
  t2 = y*y;
  t3 = z*z;
  t4 = t1+t2+t3;
  t5 = sqrt(t4);
  t7 = 1/t5/t4;
  t19 = 1/t4;
  t23 = t4*t4;
  t24 = 1/t23;

  dYdx[0] = 0.0;
  dYdx[1] = -0.488603*y*t7*x;
  dYdx[2] = -0.488603*z*t7*x;
  dYdx[3] = 0.488603/t5-0.488603*t1*t7;
  dYdx[4] = 0.1092548E1*y*t19-0.2185096E1*t1*y*t24;
  dYdx[5] = -0.2185096E1*x*y*t24*z;
  dYdx[6] = -0.1892352E1*t3*t24*x;
  dYdx[7] = 0.1092548E1*z*t19-0.2185096E1*z*t1*t24;
  dYdx[8] = 0.1092548E1*x*t19-0.1092548E1*(t1-t2)*t24*x;
}

void
GradientSH::computedYdy( Vec3f dir, float *dYdy ) {
  float t1, t2, t3, t4, t5, t9, t13, t18, t22, t23;
  float x = dir[0];
  float y = dir[1];
  float z = dir[2];

  t1 = x*x;
  t2 = y*y;
  t3 = z*z;
  t4 = t1+t2+t3;
  t5 = sqrt(t4);
  t9 = 1/t5/t4;
  t13 = y*t9;
  t18 = 1/t4;
  t22 = t4*t4;
  t23 = 1/t22;

  dYdy[0] = 0.0;
  dYdy[1] = 0.488603/t5-0.488603*t2*t9;
  dYdy[2] = -0.488603*t13*z;
  dYdy[3] = -0.488603*t13*x;
  dYdy[4] = 0.1092548E1*x*t18-0.2185096E1*x*t2*t23;
  dYdy[5] = 0.1092548E1*z*t18-0.2185096E1*t2*z*t23;
  dYdy[6] = -0.1892352E1*y*t3*t23;
  dYdy[7] = -0.2185096E1*x*y*t23*z;
  dYdy[8] = -0.1092548E1*y*t18-0.1092548E1*(t1-t2)*t23*y;
}

void
GradientSH::computedYdz( Vec3f dir, float *dYdz ) {
  float t1, t2, t3, t4, t5, t7, t20, t21, t25;
  float x = dir[0];
  float y = dir[1];
  float z = dir[2];

  t1 = x*x;
  t2 = y*y;
  t3 = z*z;
  t4 = t1+t2+t3;
  t5 = sqrt(t4);
  t7 = 1/t5/t4;
  t20 = t4*t4;
  t21 = 1/t20;
  t25 = 1/t4;

  dYdz[0] = 0.0;
  dYdz[1] = -0.488603*y*t7*z;
  dYdz[2] = 0.488603/t5-0.488603*t3*t7;
  dYdz[3] = -0.488603*z*t7*x;
  dYdz[4] = -0.2185096E1*x*y*t21*z;
  dYdz[5] = 0.1092548E1*y*t25-0.2185096E1*y*t3*t21;
  dYdz[6] = 0.1892352E1*z*t25-0.1892352E1*t3*z*t21;
  dYdz[7] = 0.1092548E1*x*t25-0.2185096E1*t3*t21*x;
  dYdz[8] = -0.1092548E1*(t1-t2)*t21*z;
}


// ------------------------------------------------
// Rendering with Gradients
// ------------------------------------------------

void 
renderSHGradLighting() {
  Vec3f position(0.f,0.f,0.f);
  SimpleSphericalHarmonics sh( *cubemapsh );	

  // do rotation of SH  (if it wasn't generated on-the-fly, in which case it's in correct coord system anyway)
  if( !generateSH ) {
	sh.rotate(Matrix4x4f(lighti.tb.getMatrix()));
	Matrix4x4f m(rotateobjecti.tb.getMatrix());
	m.transpone(); // invert, but it's only a rotation, so transpone is enough
	sh.rotate(m);
  }

  SimpleObject *obj = so;

  vector<Vec3f> colors;
  for( int vert = 0; vert < obj->getNumberVertices(); ++vert ) {
	Vec3f v = obj->getVertex( vert );
	Vec3f n = obj->getNormal( vert );

	// reset to original
	sh = *cubemapsh;

	// extrapolate according to gradient
	const int r = 0, g = 1, b = 2;
	float dx = v[0] - position[0]; //dx = -dx;
	float dy = v[1] - position[1]; //dy = -dy;
	float dz = v[2] - position[2]; //dz = -dz;
	for( int i = 0; i < (Lmax+1)*(Lmax+1); ++i ) {
	  (sh.getCoefficients())[i][r] += gradientshX.getCoefficients()[i][r]*dx + gradientshY.getCoefficients()[i][r]*dy + gradientshZ.getCoefficients()[i][r]*dz;
	  (sh.getCoefficients())[i][g] += gradientshX.getCoefficients()[i][g]*dx + gradientshY.getCoefficients()[i][g]*dy + gradientshZ.getCoefficients()[i][g]*dz;
	  (sh.getCoefficients())[i][b] += gradientshX.getCoefficients()[i][b]*dx + gradientshY.getCoefficients()[i][b]*dy + gradientshZ.getCoefficients()[i][b]*dz;
	}
	sh.scale( globalLightScale ); // some user-defined scale for the lighting

	Vec3f col;
	if( !usePRT ) {
	  // for diffuse lighting without shadowing (no PRT)
	  sh.convolveDiffuse();
	  col = sh.getRadiance( n );
	} else {
	  // using PRT
	  Vec3f *coeff = sh.getCoefficients();
	  for( int i = 9; i < 25; ++i ) coeff[i] = Vec3f(0.f,0.f,0.f); // clear higher-order stuff, we don't have derivatives...
	  col = sh.dot( occlusionFactorsSH[vert] );
	}

	colors.push_back(col);
  }

  for( int triID = 0; triID < obj->getNumberTriangles(); ++triID ) {
	int index[3];
	Vec3f v[3];

	index[0] = obj->getTriangleObject(triID)->vi[0];
	index[1] = obj->getTriangleObject(triID)->vi[1];
	index[2] = obj->getTriangleObject(triID)->vi[2];

	v[0] = obj->getVertex( index[0] );
	v[1] = obj->getVertex( index[1] );
	v[2] = obj->getVertex( index[2] );
  	
	glBegin( GL_TRIANGLES );
	glColor3fv( colors[index[0]] );
	//glTexCoord3fv( obj->getTexCoords(index[0]) );
	glTexCoord3fv( obj->getFaceTexCoords( triID, 0 ) );
	glVertex3fv( v[0] );

	glColor3fv( colors[index[1]] );
	//glTexCoord3fv( obj->getTexCoords(index[1]) );
	glTexCoord3fv( obj->getFaceTexCoords( triID, 1 ) );
	glVertex3fv( v[1] );

	glColor3fv( colors[index[2]] );
	//glTexCoord3fv( obj->getTexCoords(index[2]) );
	glTexCoord3fv( obj->getFaceTexCoords( triID, 2 ) );
	glVertex3fv( v[2] );
	glEnd();
  }
}