Got a lot of different distributions using the SG distribution

sources/plugins/nonlinear_function_spherical_gaussian/function.cpp

@@ -23,8 +23,7 @@ vec spherical_gaussian_function::operator()(const vec& x) const
 vec spherical_gaussian_function::value(const vec& x) const 
 {
 	vec res(dimY());
-	double dot;
-	params::convert(&x[0], params::CARTESIAN, params::COS_TH, &dot);
+	double dot = compute_dot(x);
 
 	for(int i=0; i<dimY(); ++i)
 	{
@@ -43,58 +42,135 @@ void spherical_gaussian_function::setDimY(int nY)
     _n.resize(_nY) ;
     _ks.resize(_nY) ;
 }
+		
+double spherical_gaussian_function::compute_dot(const vec& x) const
+{
+	double dot;
+	if(_type == Half)
+	{
+		params::convert(&x[0], params::CARTESIAN, params::COS_TH, &dot);
+	}
+	else if(_type == Back)
+	{
+		params::convert(&x[0], params::CARTESIAN, params::COS_TK, &dot);
+	}
+	else if(_type == Moment)
+	{
+		double spherical[4];
+		params::convert(&x[0], params::CARTESIAN, params::SPHERICAL_TL_PL_TV_PV, &spherical[0]);
+		spherical[1] *= _a;
+		spherical[1] += M_PI;
+		
+		double cart[6];
+		params::convert(&spherical[0], params::SPHERICAL_TL_PL_TV_PV, params::CARTESIAN, &cart[0]);
+		dot = cart[0]*cart[3] + cart[1]*cart[4] + cart[2]*cart[5];
+	}
+	else
+	{
+		double spherical[4];
+		params::convert(&x[0], params::CARTESIAN, params::SPHERICAL_TL_PL_TV_PV, &spherical[0]);
+		spherical[1] += M_PI;
+		
+		double cart[6];
+		params::convert(&spherical[0], params::SPHERICAL_TL_PL_TV_PV, params::CARTESIAN, &cart[0]);
+		dot = cart[0]*cart[3] + cart[1]*cart[4] + cart[2]*cart[5];
+	}
+
+	return dot;
+}
 
 //! Number of parameters to this non-linear function
 int spherical_gaussian_function::nbParameters() const 
 {
-    return 2*dimY();
+	if(_type == Moment)
+	{
+		return 2*dimY()+1;
+	}
+	else
+	{
+		return 2*dimY();
+	}
 }
 
 //! Get the vector of parameters for the function
 vec spherical_gaussian_function::parameters() const 
 {
-    vec res(2*dimY());
-    for(int i=0; i<dimY(); ++i)
-    {
-        res[i*2 + 0] = _ks[i];
-        res[i*2 + 1] = _n[i];
-    }
-
-    return res;
+	if(_type == Moment)
+	{
+		vec res(2*dimY()+1);
+		res[2*dimY()] = _a;
+		for(int i=0; i<dimY(); ++i)
+		{
+			res[i*2 + 0] = _ks[i];
+			res[i*2 + 1] = _n[i];
+		}
+
+		return res;
+	}
+	else
+	{
+		vec res(2*dimY());
+		for(int i=0; i<dimY(); ++i)
+		{
+			res[i*2 + 0] = _ks[i];
+			res[i*2 + 1] = _n[i];
+		}
+		return res;
+	}
 }
 
 //! \brief get the min values for the parameters
 vec spherical_gaussian_function::getParametersMin() const
 {
-    vec res(2*dimY());
-    for(int i=0; i<dimY(); ++i)
-    {
-        res[i*2 + 0] = 0.0;
-        res[i*2 + 1] = 0.0;
-    }
-
-    return res;
+	if(_type == Moment)
+	{
+		vec res(2*dimY()+1);
+		res[2*dimY()] = 0.0;
+		for(int i=0; i<dimY(); ++i)
+		{
+			res[i*2 + 0] = 0.0;
+			res[i*2 + 1] = 0.0;
+		}
+
+		return res;
+	}
+	else
+	{
+		vec res(2*dimY());
+		for(int i=0; i<dimY(); ++i)
+		{
+			res[i*2 + 0] = 0.0;
+			res[i*2 + 1] = 0.0;
+		}
+
+		return res;
+	}
 }
 
 //! Update the vector of parameters for the function
 void spherical_gaussian_function::setParameters(const vec& p) 
 {
-    for(int i=0; i<dimY(); ++i)
-    {
-        _ks[i] = p[i*2 + 0];
-        _n[i]  = p[i*2 + 1];
-    }
+	if(_type == Moment)
+	{
+		_a = p[2*dimY()];
+	}
+
+	for(int i=0; i<dimY(); ++i)
+	{
+		_ks[i] = p[i*2 + 0];
+		_n[i]  = p[i*2 + 1];
+	}
 }
 
 //! Obtain the derivatives of the function with respect to the 
 //! parameters. 
 vec spherical_gaussian_function::parametersJacobian(const vec& x) const 
 {
-	double dot;
-	params::convert(&x[0], params::CARTESIAN, params::COS_TH, &dot);
+	double dot = compute_dot(x);
 
     vec jac(dimY()*nbParameters());
 	 for(int i=0; i<dimY(); ++i)
+	 {
 		 for(int j=0; j<dimY(); ++j)
 		 {
 			 if(i == j)
@@ -105,6 +181,7 @@ vec spherical_gaussian_function::parametersJacobian(const vec& x) const
 
 				 // df / dN
 				 jac[i*nbParameters() + j*2+1] = _ks[j] * (dot-1) * exp(_n[j]*(dot - 1));
+
 			 }
 			 else
 			 {
@@ -113,6 +190,21 @@ vec spherical_gaussian_function::parametersJacobian(const vec& x) const
 			 }
 		 }
 
+		 if(_type == Moment)
+		 {
+			 // df / da
+			 double spherical[4];
+			 params::convert(&x[0], params::CARTESIAN, params::SPHERICAL_TL_PL_TV_PV, &spherical[0]);
+			 spherical[1] += 0.5* M_PI;
+
+			 double cart[6];
+			 params::convert(&spherical[0], params::SPHERICAL_TL_PL_TV_PV, params::CARTESIAN, &cart[0]);
+			 double dot2 = cart[0]*cart[3] + cart[1]*cart[4] + cart[2]*cart[5];
+			 jac[i*nbParameters()] = _ks[i] * _n[i] * exp(_n[i] * (dot - 1)) * _a * dot2;
+
+		 }
+	 }
+
     return jac;
 }
 		
@@ -123,6 +215,33 @@ void spherical_gaussian_function::bootstrap(const data* d, const arguments& args
 		_ks[i] = 1.0;
 		_n[i]  = 1.0;
 	}
+
+	// Parse the lobe type
+	if(args.is_defined("sg-type"))
+	{
+		std::string stype = args.get_string("sg-type", "half");
+		if(stype == "mirror")
+		{
+			_type = Mirror;
+		}
+		else if(stype == "half")
+		{
+			_type = Half;
+		}
+		else if(stype == "back")
+		{
+			_type = Back;
+		}
+		else if(stype == "retro")
+		{
+			_type = Retro;
+		}
+		else if(stype == "moment")
+		{
+			_type = Moment;
+			_a = 1.0;
+		}
+	}
 }
 
 //! Load function specific files
@@ -156,6 +275,11 @@ void spherical_gaussian_function::load(std::istream& in)
 		in >> token >> _ks[i];
 		in >> token >> _n[i];
 	}
+	
+	if(_type == Moment)
+	{
+		in >> token >> _a;
+	}
 }
 
 
@@ -172,6 +296,11 @@ void spherical_gaussian_function::save_call(std::ostream& out, const arguments&
 			 out << "Ks " << _ks[i] << std::endl;
 			 out << "N  " <<  _n[i] << std::endl;
 		 }
+	
+		 if(_type == Moment)
+		 {
+			 out << "a  " << _a;
+		 }
 
 		 out << std::endl;
 	 }
@@ -204,8 +333,37 @@ void spherical_gaussian_function::save_body(std::ostream& out, const arguments&
     {
         out << "vec3 spherical_gaussian(vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y, vec3 ks, vec3 Nl)" << std::endl;
         out << "{" << std::endl;
-		  out << "\tvec3 H = normalize(L + V);" << std::endl;
-        out << "\tvec3 ext_dot = vec3(dot(H,N));" << std::endl;
+
+		  switch(_type)
+		  {
+			  case Half:
+				  out << "\tvec3 H = normalize(L + V);" << std::endl;
+				  out << "\tvec3 ext_dot = vec3(dot(H,N));" << std::endl;
+				  break;
+
+			  case Moment:
+				  out << "\tfloat t = " << _a << "*acos(dot(V, N));" << std::endl;
+				  out << "\tfloat p = atan(dot(V, Y), dot(V, X)) + " << M_PI << ";" << std::endl;
+				  out << "\tvec3 R = cos(p)*sin(t)*X + sin(p)*sin(t)*Y + cos(t)*N;" << std::endl;
+				  out << "\tvec3 ext_dot = vec3(dot(L,R));" << std::endl;
+				  break;
+
+			  case Back:
+				  out << "\tvec3 Vp = 2.0f*N*(dot(N,V)) - V;" << std::endl;
+				  out << "\tvec3 K  = normalize(L + Vp);" << std::endl;
+				  out << "\tvec3 ext_dot = vec3(dot(K,N));" << std::endl;
+				  break;
+			  
+			  case Retro:
+				  out << "\tvec3 ext_dot = vec3(dot(L,V));" << std::endl;
+				  break;
+
+			  default:
+				  out << "\tvec3 R = 2*dot(V,N)*N - V;" << std::endl;
+				  out << "\tvec3 ext_dot = vec3(dot(L,R));" << std::endl;
+				  break;
+		  }
+
         out << "\treturn ks * exp(Nl * (ext_dot - vec3(1)));" << std::endl;
         out << "}" << std::endl;
     }

sources/plugins/nonlinear_function_spherical_gaussian/function.h

@@ -31,7 +31,16 @@ class spherical_gaussian_function : public nonlinear_function
 
 	public: // methods
 
-		spherical_gaussian_function() : _n(1) 
+		enum type
+		{
+			Mirror,
+			Half,
+			Retro,
+			Back,
+			Moment
+		};
+
+		spherical_gaussian_function() : _a(1), _type(Mirror)
 		{ 
 			setParametrization(params::CARTESIAN);
 			setDimX(6);
@@ -78,14 +87,18 @@ class spherical_gaussian_function : public nonlinear_function
 		//! \brief Set the number of output dimensions
 		void setDimY(int nY);
 
-		//! \brief Set the number of lobes to be used in the fit
-		void setNbLobes(int N);
-
 	private: // methods
 
+		//! \brief Compute the cosine for inside the lobe function. 
+		//! Depending on the lobe type, the dot product can have 
+		//! different evaluations.
+		double compute_dot(const vec& in) const;
 
 	private: // data
 
 		vec _n, _ks; // Lobes data
+		double _a;   // Scaling factor for the moment based SG
+
+		type _type;  // Lobe type
 } ;