From 51dddaccacef51e6a7fb2d086ee2f62cf693fab3 Mon Sep 17 00:00:00 2001
From: Laurent Belcour <laurent.belcour@inria.fr>
Date: Tue, 22 Oct 2013 14:36:24 +0200
Subject: [PATCH] Fitting with fresnel is now working properly

---
 sources/core/args.h                           |   2 +-
 sources/core/common.cpp                       |   4 +-
 sources/core/function.cpp                     | 371 ++++++++++++++++--
 sources/core/function.h                       | 271 +------------
 sources/core/vertical_segment.cpp             |   2 +-
 .../plugins/nonlinear_fitter_ceres/fitter.cpp |   3 +
 .../function.cpp                              |  41 +-
 .../function.h                                |   3 +-
 .../nonlinear_fresnel_schlick/function.cpp    |   4 +
 .../nonlinear_fresnel_schlick/function.h      |   1 +
 .../nonlinear_function_diffuse/function.h     |   1 +
 sources/softs/tests/main.cpp                  |  52 ++-
 12 files changed, 452 insertions(+), 303 deletions(-)

diff --git a/sources/core/args.h b/sources/core/args.h
index bd78416..5ce0949 100644
--- a/sources/core/args.h
+++ b/sources/core/args.h
@@ -307,7 +307,7 @@ class arguments
 			  }
 
 			  arguments current_args(argc, argv);
-			  delete argv;
+			  delete[] argv;
 			  return current_args;
 		  }
 
diff --git a/sources/core/common.cpp b/sources/core/common.cpp
index 3bf9578..dd013de 100644
--- a/sources/core/common.cpp
+++ b/sources/core/common.cpp
@@ -53,7 +53,7 @@ double dot(const vec& a, const vec& b)
 
 vec product(const vec& a, const vec& b)
 {
-    if(a.size() == 1 && b.size() > 1)
+    if(a.size() == 1 && b.size() >= 1)
     {
         vec res(b.size());
         for(int i=0; i<b.size(); ++i)
@@ -62,7 +62,7 @@ vec product(const vec& a, const vec& b)
         }
         return res;
     }
-    else if(b.size() == 1 && a.size() > 1)
+    else if(b.size() == 1 && a.size() >= 1)
     {
         vec res(a.size());
         for(int i=0; i<a.size(); ++i)
diff --git a/sources/core/function.cpp b/sources/core/function.cpp
index 406da25..270325a 100644
--- a/sources/core/function.cpp
+++ b/sources/core/function.cpp
@@ -308,22 +308,50 @@ vec compound_function::operator()(const vec& x) const
 }
 vec compound_function::value(const vec& x) const
 {
-	vec res(_nY);
-	res = vec::Zero(_nY);
+	vec res = vec::Zero(dimY());
 	for(unsigned int i=0; i<fs.size(); ++i)
 	{
-		if(fs[i]->input_parametrization() != input_parametrization())
-		{
-			vec temp_x(fs[i]->dimX());
-			params::convert(&x[0], input_parametrization(), fs[i]->input_parametrization(), &temp_x[0]);
-			res = res + fs[i]->value(temp_x);
-		}
-		else
+		vec temp_x(fs[i]->dimX());
+		params::convert(&x[0], input_parametrization(), fs[i]->input_parametrization(), &temp_x[0]);
+		res = res + fs[i]->value(temp_x);
+	}
+	return res;
+}
+
+vec compound_function::parametersJacobian(const vec& x) const
+{
+	int nb_params = nbParameters();
+	vec jac = vec::Zero(nb_params*dimY());
+
+	int start_i = 0;
+
+	// Export the sub-Jacobian for each function
+	for(unsigned int f=0; f<fs.size(); ++f)
+	{
+		nonlinear_function* func = fs[f];
+		int nb_f_params = func->nbParameters(); 
+
+		// Only export Jacobian if there are non-linear parameters
+		if(nb_f_params > 0 && !is_fixed[f])
 		{
-			res = res + fs[i]->value(x);
+
+			vec temp_x(func->dimX());
+			params::convert(&x[0], input_parametrization(), func->input_parametrization(), &temp_x[0]);
+			vec func_jac = func->parametersJacobian(temp_x);
+
+			for(int i=0; i<nb_f_params; ++i)
+			{
+				for(int y=0; y<dimY(); ++y)
+				{
+					jac[y*nb_params + (i+start_i)] = func_jac[y*nb_f_params + i];
+				}
+			}
+
+			start_i += nb_f_params;
 		}
 	}
-	return res;
+
+	return jac;
 }
 
 void compound_function::push_back(nonlinear_function* f, const arguments& f_args)
@@ -361,6 +389,62 @@ nonlinear_function* compound_function::operator[](int i) const
 #endif
 	return fs[i];
 }
+		
+unsigned int compound_function::size() const
+{
+	return fs.size();
+}
+
+bool compound_function::load(std::istream& in)
+{
+	int nb_good = 0; // Number of correctly openned functions
+	for(unsigned int i=0; i<fs.size(); ++i)
+	{
+		std::streampos pos = in.tellg();
+		if(!fs[i]->load(in))
+		{
+			in.seekg(pos);
+		}
+		else
+		{
+			nb_good++;
+		}
+	}
+	return nb_good > 0;
+}
+
+void compound_function::setParametrization(params::input new_param)
+{
+	if(new_param == params::UNKNOWN_INPUT)
+		return;
+
+	// If there is more than one parametrization defined to this function, convert
+	// it to the parametrization that conserves the most degrees of freedom. Right now
+	// this is the CARTESIAN param. It might change in future release.
+	if(input_parametrization() != new_param)
+	{
+		function::setParametrization(params::CARTESIAN);
+		function::setDimX(6);
+	}
+
+	for(unsigned int i=0; i<fs.size(); ++i)
+	{
+		if(fs[i]->input_parametrization() == params::UNKNOWN_INPUT)
+		{
+			fs[i]->setParametrization(new_param);
+			fs[i]->setDimX(params::dimension(new_param));
+		}
+	}
+}
+		
+void compound_function::setParametrization(params::output new_param)
+{
+	parametrized::setParametrization(new_param);
+	for(unsigned int i=0; i<fs.size(); ++i)
+	{
+		fs[i]->setParametrization(new_param);
+	}
+}
 
 void compound_function::bootstrap(const ::data* d, const arguments& args)
 {
@@ -404,6 +488,11 @@ void compound_function::bootstrap(const ::data* d, const arguments& args)
 				for(unsigned int i=0; i<fs.size(); ++i)
 				{
 					std::streampos pos = file.tellg();
+					/*
+					if(dynamic_cast<product_function*>(fs[i]) != NULL)
+					{
+					}
+					*/
 
 					// If the function cannot be loaded, put the input stream
 					// in the previous state and bootstrap normaly this function.
@@ -431,6 +520,200 @@ void compound_function::bootstrap(const ::data* d, const arguments& args)
 	}
 }
 
+void compound_function::setDimX(int nX) 
+{
+	if(input_parametrization() == params::UNKNOWN_INPUT)
+	{
+		function::setDimX(nX);
+	}
+	
+	for(unsigned int i=0; i<fs.size(); ++i)
+	{
+		if(fs[i]->input_parametrization() == params::UNKNOWN_INPUT)
+		{
+			fs[i]->setDimX(nX);
+		}
+	}
+}
+		
+void compound_function::setDimY(int nY)
+{
+	function::setDimY(nY);
+	for(unsigned int i=0; i<fs.size(); ++i)
+	{
+		fs[i]->setDimY(nY);
+	}
+}
+
+void compound_function::setMin(const vec& min) 
+{
+	function::setMin(min);
+	for(unsigned int i=0; i<fs.size(); ++i)
+	{
+		fs[i]->setMin(min);
+	}
+}
+void compound_function::setMax(const vec& max) 
+{
+	function::setMax(max);
+	for(unsigned int i=0; i<fs.size(); ++i)
+	{
+		fs[i]->setMax(max);
+	}
+}
+
+//! Number of parameters to this non-linear function
+int compound_function::nbParameters() const
+{
+	int nb_params = 0;
+	for(unsigned int i=0; i<fs.size(); ++i)
+	{
+		if(!is_fixed[i]) {
+			nb_params += fs[i]->nbParameters();
+		}
+	}
+	return nb_params;
+}
+
+//! Get the vector of parameters for the function
+vec compound_function::parameters() const
+{
+	vec params(nbParameters());
+	int current_i = 0;
+	for(unsigned int f=0; f<fs.size(); ++f)
+	{
+		int f_size = fs[f]->nbParameters();
+
+		// Handle when there is no parameters to include
+		if(f_size > 0 && !is_fixed[f])
+		{
+			vec f_params = fs[f]->parameters();
+			for(int i=0; i<f_size; ++i)
+			{
+				params[i + current_i] = f_params[i];
+			}
+
+			current_i += f_size;
+		}
+	}
+
+	return params;
+}
+
+//! Get the vector of min parameters for the function
+vec compound_function::getParametersMin() const
+{
+	vec params(nbParameters());
+	int current_i = 0;
+	for(unsigned int f=0; f<fs.size(); ++f)
+	{
+		int f_size = fs[f]->nbParameters();
+
+		// Handle when there is no parameters to include
+		if(f_size > 0 && !is_fixed[f])
+		{
+			vec f_params = fs[f]->getParametersMin();
+			for(int i=0; i<f_size; ++i)
+			{
+				params[i + current_i] = f_params[i];
+			}
+
+			current_i += f_size;
+		}
+	}
+
+	return params;
+}
+
+//! Get the vector of min parameters for the function
+vec compound_function::getParametersMax() const
+{
+	vec params(nbParameters());
+	int current_i = 0;
+	for(unsigned int f=0; f<fs.size(); ++f)
+	{
+		int f_size = fs[f]->nbParameters();
+
+		// Handle when there is no parameters to include
+		if(f_size > 0 && !is_fixed[f])
+		{
+			vec f_params = fs[f]->getParametersMax();
+			for(int i=0; i<f_size; ++i)
+			{
+				params[i + current_i] = f_params[i];
+			}
+
+			current_i += f_size;
+		}
+	}
+
+	return params;
+}
+
+//! Update the vector of parameters for the function
+void compound_function::setParameters(const vec& p) 
+{
+	int current_i = 0;
+	for(unsigned int f=0; f<fs.size(); ++f)
+	{
+		int f_size = fs[f]->nbParameters();
+
+		// Handle when there is no parameters to include
+		if(f_size > 0 && !is_fixed[f])
+		{
+			vec f_params(f_size);
+			for(int i=0; i<f_params.size(); ++i)
+			{
+				f_params[i] = p[i + current_i];
+			}
+
+			fs[f]->setParameters(f_params);
+			current_i += f_size;
+		}
+	}
+}
+		
+void compound_function::save_body(std::ostream& out, const arguments& args) const
+{
+	for(unsigned int i=0; i<fs.size(); ++i)
+	{
+		fs[i]->save_body(out, args);
+		out << std::endl;
+	}
+
+	function::save_body(out, args);
+}
+		
+void compound_function::save_call(std::ostream& out, const arguments& args) const
+{
+	bool is_cpp    = args["export"] == "C++";
+	bool is_shader = args["export"] == "shader" || args["export"] == "explorer";
+	bool is_matlab = args["export"] == "matlab";
+
+	// This part is export specific. For ALTA, the coefficients are just
+	// dumped as is with a #FUNC {plugin_name} header.
+	//
+	// For C++ export, the function call should be done before hand and
+	// the line should look like:
+	//   res += call_i(x);
+	for(unsigned int i=0; i<fs.size(); ++i)
+	{
+		if(i != 0 && (is_cpp || is_matlab || is_shader))
+		{
+			out << "\tres += ";
+		}
+
+		fs[i]->save_call(out, args);
+
+		if(is_cpp || is_matlab || is_shader)
+		{
+			out << ";" << std::endl;
+		}
+	}
+
+	function::save_call(out, args);
+}
+
 
 
 /*--- Product functions implementation ----*/
@@ -439,6 +722,19 @@ product_function::product_function(nonlinear_function* g1, nonlinear_function* g
 {
 	this->f1 = g1;
 	this->f2 = g2;
+
+	// If the two parametrization are different, use the CARTESIAN parametrization
+	// as the input parametrization, then do the convertion for all the functions.
+	if(g1->input_parametrization() != g2->input_parametrization())
+	{
+		function::setParametrization(params::CARTESIAN);
+		function::setDimX(6);
+	}
+	else
+	{
+		setParametrization(g1->input_parametrization());
+		function::setDimX(g1->dimX());
+	}
 }
 
 vec product_function::operator()(const vec& x) const
@@ -447,15 +743,23 @@ vec product_function::operator()(const vec& x) const
 }
 vec product_function::value(const vec& x) const
 {
-	// Get the first function value
-	vec f1res = f1->value(x);
-
-	// Convert input space to the 2nd function parametreization
+	// Convert input space to the 1rst function parametrization and compute the
+	// output value
+	vec xf1(f1->dimX());
+	params::convert(&x[0], input_parametrization(), f1->input_parametrization(), &xf1[0]);
+	vec f1res = f1->value(xf1);
+
+	// Convert input space to the 2nd function parametrization and compute the
+	// output value
 	vec xf2(f2->dimX());
-	params::convert(&x[0], f1->input_parametrization(), f2->input_parametrization(), &xf2[0]);
+	params::convert(&x[0], input_parametrization(), f2->input_parametrization(), &xf2[0]);
 	vec f2res = f2->value(xf2);
 
-	return product(f1res, f2res);
+	vec res = product(f1res, f2res);
+
+	//std::cout << f1res << " :::::::::: " << res << std::endl;
+	//std::cout << f2res << std::endl;
+	return res;
 }
 		
 bool product_function::load(std::istream& in)
@@ -568,6 +872,7 @@ void product_function::bootstrap(const data* d, const arguments& args)
 
 					// Bootstrap the function as if it was not loaded
 					f1->bootstrap(d, args);
+					std::cout << "<<DEBUG>> Unable to load first function of product, regular bootstraping" << std::endl;
 				}
 				
 				// If the second function cannot be loaded, put the input stream
@@ -578,6 +883,7 @@ void product_function::bootstrap(const data* d, const arguments& args)
 
 					// Bootstrap the function as if it was not loaded
 					f2->bootstrap(d, args);
+					std::cout << "<<DEBUG>> Unable to load second function of product, regular bootstraping" << std::endl;
 				}
 			}
 			else
@@ -726,14 +1032,16 @@ vec product_function::parametersJacobian(const vec& x) const
 	int nb_params = nb_f1_params + nb_f2_params;
 
 	// Convert the input value x to the input space of the f1tion
-	vec xf(f2->dimX());
-	params::convert(&x[0], f1->input_parametrization(), f2->input_parametrization(), &xf[0]);
+	vec xf2(f2->dimX());
+	params::convert(&x[0], input_parametrization(), f2->input_parametrization(), &xf2[0]);
+	vec xf1(f1->dimX());
+	params::convert(&x[0], input_parametrization(), f1->input_parametrization(), &xf1[0]);
 
-	vec f1_jacobian = f1->parametersJacobian(x);
-	vec f2_jacobian = f2->parametersJacobian(xf);
+	vec f1_jacobian = f1->parametersJacobian(xf1);
+	vec f2_jacobian = f2->parametersJacobian(xf2);
 
-	vec f1_value = f1->value(x);
-	vec f2_value = f2->value(xf);
+	vec f1_value = f1->value(xf1);
+	vec f2_value = f2->value(xf2);
 
 	// F = f2nel; f = f1tion
 	// d(F * f)(x) /dp = F(x) df(x) /dp + f(x) dF(x) / dp
@@ -753,11 +1061,6 @@ vec product_function::parametersJacobian(const vec& x) const
 
 	return jac;
 }
-		
-params::input product_function::input_parametrization() const
-{
-	return f1->input_parametrization();
-}
 
 //! \brief provide the outout parametrization of the object.
 params::output product_function::output_parametrization() const
@@ -767,8 +1070,16 @@ params::output product_function::output_parametrization() const
 		
 void product_function::setParametrization(params::input new_param)
 {
-	f1->setParametrization(new_param);
-	f2->setParametrization(new_param);
+	if(f1->input_parametrization() == params::UNKNOWN_INPUT) 
+	{
+		f1->setParametrization(new_param); 
+		f1->setDimX(params::dimension(new_param));
+	}
+	if(f2->input_parametrization() == params::UNKNOWN_INPUT) 
+	{
+		f2->setParametrization(new_param);
+		f2->setDimX(params::dimension(new_param));
+	}
 }
 
 void product_function::setParametrization(params::output new_param)
diff --git a/sources/core/function.h b/sources/core/function.h
index f5b840b..289493e 100644
--- a/sources/core/function.h
+++ b/sources/core/function.h
@@ -201,29 +201,10 @@ class compound_function: public nonlinear_function
 		nonlinear_function* operator[](int i) const;
 
 		//! \brief Access to the number of elements in the compound object.
-		unsigned int size() const
-		{
-			return fs.size();
-		}
+		unsigned int size() const;
 
 		//! Load function specific files
-		virtual bool load(std::istream& in)
-		{
-			int nb_good = 0; // Number of correctly openned functions
-			for(unsigned int i=0; i<fs.size(); ++i)
-			{
-				std::streampos pos = in.tellg();
-				if(!fs[i]->load(in))
-				{
-					in.seekg(pos);
-				}
-				else
-				{
-					nb_good++;
-				}
-			}
-			return nb_good > 0;
-		}
+		virtual bool load(std::istream& in);
 
 		//! \brief Provide a first rough fit of the function. 
 		//! For compound object, you can define the first guess using the
@@ -248,152 +229,29 @@ class compound_function: public nonlinear_function
 		virtual void bootstrap(const ::data* d, const arguments& args);
 
 		//! Set the dimension of the input space of the function
-		virtual void setDimX(int nX) 
-		{
-			function::setDimX(nX);
-			for(unsigned int i=0; i<fs.size(); ++i)
-			{
-				fs[i]->setDimX(nX);
-			}
-		}
+		virtual void setDimX(int nX);
+
 		//! Set the dimension of the output space of the function
-		virtual void setDimY(int nY)
-		{
-			function::setDimY(nY);
-			for(unsigned int i=0; i<fs.size(); ++i)
-			{
-				fs[i]->setDimY(nY);
-			}
-		}
+		virtual void setDimY(int nY);
 
 		// Acces to the domain of definition of the function
-		virtual void setMin(const vec& min) 
-		{
-			function::setMin(min);
-			for(unsigned int i=0; i<fs.size(); ++i)
-			{
-				fs[i]->setMin(min);
-			}
-		}
-		virtual void setMax(const vec& max) 
-		{
-			function::setMax(max);
-			for(unsigned int i=0; i<fs.size(); ++i)
-			{
-				fs[i]->setMax(max);
-			}
-		}
+		virtual void setMin(const vec& min);
+		virtual void setMax(const vec& max);
 
 		//! Number of parameters to this non-linear function
-		virtual int nbParameters() const
-		{
-			int nb_params = 0;
-			for(unsigned int i=0; i<fs.size(); ++i)
-			{
-				if(!is_fixed[i]) {
-					nb_params += fs[i]->nbParameters();
-				}
-			}
-			return nb_params;
-		}
+		virtual int nbParameters() const;
 
 		//! Get the vector of parameters for the function
-		virtual vec parameters() const
-		{
-			vec params(nbParameters());
-			int current_i = 0;
-			for(unsigned int f=0; f<fs.size(); ++f)
-			{
-				int f_size = fs[f]->nbParameters();
-
-				// Handle when there is no parameters to include
-				if(f_size > 0 && !is_fixed[f])
-				{
-					vec f_params = fs[f]->parameters();
-					for(int i=0; i<f_size; ++i)
-					{
-						params[i + current_i] = f_params[i];
-					}
-
-					current_i += f_size;
-				}
-			}
-
-			return params;
-		}
+		virtual vec parameters() const;
 		
 		//! Get the vector of min parameters for the function
-		virtual vec getParametersMin() const
-		{
-			vec params(nbParameters());
-			int current_i = 0;
-			for(unsigned int f=0; f<fs.size(); ++f)
-			{
-				int f_size = fs[f]->nbParameters();
-
-				// Handle when there is no parameters to include
-				if(f_size > 0 && !is_fixed[f])
-				{
-					vec f_params = fs[f]->getParametersMin();
-					for(int i=0; i<f_size; ++i)
-					{
-						params[i + current_i] = f_params[i];
-					}
-
-					current_i += f_size;
-				}
-			}
-
-			return params;
-		}
+		virtual vec getParametersMin() const;
 		
 		//! Get the vector of min parameters for the function
-		virtual vec getParametersMax() const
-		{
-			vec params(nbParameters());
-			int current_i = 0;
-			for(unsigned int f=0; f<fs.size(); ++f)
-			{
-				int f_size = fs[f]->nbParameters();
-
-				// Handle when there is no parameters to include
-				if(f_size > 0 && !is_fixed[f])
-				{
-					vec f_params = fs[f]->getParametersMax();
-					for(int i=0; i<f_size; ++i)
-					{
-						params[i + current_i] = f_params[i];
-					}
-
-					current_i += f_size;
-				}
-			}
-
-			return params;
-		}
+		virtual vec getParametersMax() const;
 
 		//! Update the vector of parameters for the function
-		virtual void setParameters(const vec& p) 
-		{
-			int current_i = 0;
-			for(unsigned int f=0; f<fs.size(); ++f)
-			{
-				int f_size = fs[f]->nbParameters();
-
-				// Handle when there is no parameters to include
-				if(f_size > 0 && !is_fixed[f])
-				{
-					vec f_params(f_size);
-					for(int i=0; i<f_params.size(); ++i)
-					{
-						f_params[i] = p[i + current_i];
-					}
-
-					fs[f]->setParameters(f_params);
-					current_i += f_size;
-				}
-			}
-		}
+		virtual void setParameters(const vec& p);
 
 		//! \brief Obtain the derivatives of the function with respect to the 
 		//! parameters. 
@@ -404,118 +262,26 @@ class compound_function: public nonlinear_function
 		//
 		// The result vector should be orderer as res[i + dimY()*j], output
 		// dimension first, then parameters.
-		virtual vec parametersJacobian(const vec& x) const
-		{
-			int nb_params = nbParameters();
-			vec jac(nb_params*_nY);
-			jac = vec::Zero(nb_params*_nY);
-
-			int start_i = 0;
-
-			// Export the sub-Jacobian for each function
-			for(unsigned int f=0; f<fs.size(); ++f)
-			{
-				nonlinear_function* func = fs[f];
-				int nb_f_params = func->nbParameters(); 
-
-				// Only export Jacobian if there are non-linear parameters
-				if(nb_f_params > 0 && !is_fixed[f])
-				{
-
-					vec temp_x(func->dimX());
-					params::convert(&x[0], input_parametrization(), func->input_parametrization(), &temp_x[0]);
-					vec func_jac = func->parametersJacobian(temp_x);
-
-					for(int i=0; i<nb_f_params; ++i)
-					{
-						for(int y=0; y<_nY; ++y)
-						{
-                            jac[y*nb_params + (i+start_i)] = func_jac[y*nb_f_params + i];
-						}
-					}
-
-					start_i += nb_f_params;
-				}
-			}
-
-			return jac;
-		}
+		virtual vec parametersJacobian(const vec& x) const;
 
 		//! \brief can set the input parametrization of a non-parametrized
 		//! object. Print an error if it is already defined.
-		virtual void setParametrization(params::input new_param)
-		{
-			if(new_param == params::UNKNOWN_INPUT)
-				return;
-
-			parametrized::setParametrization(new_param);
-			for(unsigned int i=0; i<fs.size(); ++i)
-			{
-				if(fs[i]->input_parametrization() == params::UNKNOWN_INPUT)
-				{
-					fs[i]->setParametrization(new_param);
-				}
-			}
-		}
+		virtual void setParametrization(params::input new_param);
 
 		//! \brief can set the output parametrization of a non-parametrized
 		//! function. Throw an exception if it tries to erase a previously
 		//! defined one.
-		virtual void setParametrization(params::output new_param)
-		{
-			parametrized::setParametrization(new_param);
-			for(unsigned int i=0; i<fs.size(); ++i)
-			{
-				fs[i]->setParametrization(new_param);
-			}
-		}
+		virtual void setParametrization(params::output new_param);
 
 		//! \brief save function specific data. This has no use for ALTA export
 		//! but allows to factorize the code in the C++ or matlab export by
 		//! defining function calls that are common to all the plugins.
-		virtual void save_body(std::ostream& out, const arguments& args) const
-		{
-			for(unsigned int i=0; i<fs.size(); ++i)
-			{
-				fs[i]->save_body(out, args);
-				out << std::endl;
-			}
-
-			function::save_body(out, args);
-		}
+		virtual void save_body(std::ostream& out, const arguments& args) const;
 
 		//! \brief save object specific information. For an ALTA export the
 		//! coefficients will be exported. For a C++ or matlab export, the call
 		//! to the associated function will be done.
-		virtual void save_call(std::ostream& out, const arguments& args) const
-		{
-			bool is_cpp    = args["export"] == "C++";
-			bool is_shader = args["export"] == "shader" || args["export"] == "explorer";
-			bool is_matlab = args["export"] == "matlab";
-
-			// This part is export specific. For ALTA, the coefficients are just
-			// dumped as is with a #FUNC {plugin_name} header.
-			//
-			// For C++ export, the function call should be done before hand and
-			// the line should look like:
-			//   res += call_i(x);
-			for(unsigned int i=0; i<fs.size(); ++i)
-			{
-				if(i != 0 && (is_cpp || is_matlab || is_shader))
-				{
-					out << "\tres += ";
-				}
-
-				fs[i]->save_call(out, args);
-
-				if(is_cpp || is_matlab || is_shader)
-				{
-					out << ";" << std::endl;
-				}
-			}
-
-			function::save_call(out, args);
-		}
+		virtual void save_call(std::ostream& out, const arguments& args) const;
 
 	protected:
 		std::vector<nonlinear_function*> fs;
@@ -578,8 +344,7 @@ class product_function : public nonlinear_function
 		virtual void setMin(const vec& min);
 		virtual void setMax(const vec& max);
 		
-		//! Provide the input/output parametrization of the object.
-		virtual params::input  input_parametrization() const;
+		//! Provide the output parametrization of the object.
 		virtual params::output output_parametrization() const;
 		
 		//! Set the input/output parametrization of a non-parametrized
diff --git a/sources/core/vertical_segment.cpp b/sources/core/vertical_segment.cpp
index ae74962..c28c294 100644
--- a/sources/core/vertical_segment.cpp
+++ b/sources/core/vertical_segment.cpp
@@ -93,7 +93,7 @@ void vertical_segment::load(const std::string& filename, const arguments& args)
         else
         {
 
-		vec v(dimX() + 3*dimY()) ;
+		vec v = vec::Zero(dimX() + 3*dimY()) ;
 		for(int i=0; i<dimX(); ++i)
 			linestream >> v[i] ;
 
diff --git a/sources/plugins/nonlinear_fitter_ceres/fitter.cpp b/sources/plugins/nonlinear_fitter_ceres/fitter.cpp
index d7371bd..ad3eb33 100644
--- a/sources/plugins/nonlinear_fitter_ceres/fitter.cpp
+++ b/sources/plugins/nonlinear_fitter_ceres/fitter.cpp
@@ -211,6 +211,9 @@ bool nonlinear_fitter_ceres::fit_data(const data* d, function* fit, const argume
          /* the following starting values provide a rough fit. */
          vec p = nf->parameters();
 
+			std::cout << "<<DEBUG>> Starting vector: " << p << std::endl;
+			std::cout << "<<DEBUG>> Final vector should be between " << nf->getParametersMin() << " and " << nf->getParametersMax() << std::endl;
+
          // Create the problem
          ceres::Problem problem;
          for(int i=0; i<d->size(); ++i)
diff --git a/sources/plugins/nonlinear_fresnel_normalized_schlick/function.cpp b/sources/plugins/nonlinear_fresnel_normalized_schlick/function.cpp
index c39fdc6..a6836ad 100644
--- a/sources/plugins/nonlinear_fresnel_normalized_schlick/function.cpp
+++ b/sources/plugins/nonlinear_fresnel_normalized_schlick/function.cpp
@@ -93,7 +93,7 @@ void schlick::save_body(std::ostream& out, const arguments& args) const
         out << "vec3 normalized_schlick_fresnel(vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y, vec3 R)" << std::endl;
 		out << "{" << std::endl;
         out << "\tvec3 H = normalize(L + V);" << std::endl;
-        out << "\treturn vec3(1.0f) + (R - vec3(1.0f)) * pow(1.0f - clamp(dot(H,V), 0.0f, 1.0f), 5);" << std::endl;
+        out << "\treturn vec3(1.0f) + (vec3(1.0f)/R - vec3(1.0f)) * pow(1.0f - clamp(dot(H,V), 0.0f, 1.0f), 5);" << std::endl;
 		out << "}" << std::endl;
         out << std::endl;
 	}
@@ -102,19 +102,25 @@ void schlick::save_body(std::ostream& out, const arguments& args) const
 
 vec schlick::value(const vec& x) const
 {
-    double xp[3], cart[6];
-    params::convert(&x[0], input_parametrization(), params::RUSIN_VH, xp);
-    params::convert(&x[0], input_parametrization(), params::CARTESIAN, cart);
+    vec xp(3), cart(6);
+    params::convert(&x[0], input_parametrization(), params::RUSIN_VH, &xp[0]);
+    params::convert(&x[0], input_parametrization(), params::CARTESIAN, &cart[0]);
+	 //std::cout << xp << " ===== " << cart << std::endl;
+
 
     const double dotVH = xp[0]*cart[0] + xp[1]*cart[1] + xp[2]*cart[2];
 
-    vec res(_nY);
-    for(int i=0; i<_nY; ++i)
+    vec res(dimY());
+    for(int i=0; i<dimY(); ++i)
     {
-        res[i] = 1.0 + (R[i] - 1.0) * pow(1.0 - clamp(dotVH, 0.0, 1.0), 5.0);
+        res[i] = 1.0 + (1.0/R[i] - 1.0) * pow(1.0 - clamp(dotVH, 0.0, 1.0), 5.0);
     }
-
-    return res;
+/*
+	 std::cout << " # " << R << ", " << dotVH << std::endl;
+	 std::cout << " ## " << res << std::endl;
+    
+	 */
+	 return res;
 }
 
 //! \brief Number of parameters to this non-linear function
@@ -123,10 +129,17 @@ int schlick::nbParameters() const
     return dimY();
 }
 
-vec schlick::parametersMin() const
+vec schlick::getParametersMin() const
+{
+    vec m(nbParameters());
+    for(int i=0; i<nbParameters(); ++i) { m[i] = 1.0E-10; }
+    return m;
+}
+
+vec schlick::getParametersMax() const
 {
-    vec m(dimY());
-    for(int i=0; i<dimY(); ++i) { m[i] = 1.0; }
+    vec m(nbParameters());
+    for(int i=0; i<nbParameters(); ++i) { m[i] = 1.0; }
     return m;
 }
 
@@ -161,7 +174,7 @@ vec schlick::parametersJacobian(const vec& x) const
         for(int j=0; j<nY; ++j)
         {
             if(i == j) {
-                jac[j*dimY() + i] = pow(1.0 - clamp(dotVH, 0.0, 1.0), 5.0);
+                jac[j*dimY() + i] = - (1.0/R[i]) * pow(1.0 - clamp(dotVH, 0.0, 1.0), 5.0);
             } else {
                 jac[j*dimY() + i] = 0.0;
             }
@@ -174,5 +187,5 @@ vec schlick::parametersJacobian(const vec& x) const
 
 void schlick::bootstrap(const data* d, const arguments& args)
 {
-    for(int i=0; i<dimY(); ++i) { R[i] = 1.0; }
+    for(int i=0; i<dimY(); ++i) { R[i] = 0.000001; }
 }
diff --git a/sources/plugins/nonlinear_fresnel_normalized_schlick/function.h b/sources/plugins/nonlinear_fresnel_normalized_schlick/function.h
index 3523115..c8d2fa7 100644
--- a/sources/plugins/nonlinear_fresnel_normalized_schlick/function.h
+++ b/sources/plugins/nonlinear_fresnel_normalized_schlick/function.h
@@ -41,7 +41,8 @@ class schlick : public nonlinear_function
 		virtual void setParameters(const vec& p) ;
 
 		//! Get the vector of min parameters for the function
-		virtual vec parametersMin() const;
+		virtual vec getParametersMin() const;
+		virtual vec getParametersMax() const;
 
 		//! \brief Obtain the derivatives of the function with respect to the
 		//! parameters. 
diff --git a/sources/plugins/nonlinear_fresnel_schlick/function.cpp b/sources/plugins/nonlinear_fresnel_schlick/function.cpp
index 01a42e1..844d85f 100644
--- a/sources/plugins/nonlinear_fresnel_schlick/function.cpp
+++ b/sources/plugins/nonlinear_fresnel_schlick/function.cpp
@@ -105,6 +105,10 @@ vec schlick::value(const vec& x) const
 
     const double dotVH = xp[0]*cart[0] + xp[1]*cart[1] + xp[2]*cart[2];
 
+	 //std::cout << "v_h  = [" << xp[0] << ", " << xp[1] << ", " << xp[2] << "]" << std::endl;
+	 //std::cout << "cart = [" << cart[0] << ", " << cart[1] << ", " << cart[2] << "], [" << cart[3] << ", " << cart[4] << ", " << cart[5] << "]" << std::endl;
+	 //std::cout << dotVH << std::endl;
+
     vec res(dimY());
     for(int i=0; i<dimY(); ++i)
     {
diff --git a/sources/plugins/nonlinear_fresnel_schlick/function.h b/sources/plugins/nonlinear_fresnel_schlick/function.h
index a628cd4..cffea78 100644
--- a/sources/plugins/nonlinear_fresnel_schlick/function.h
+++ b/sources/plugins/nonlinear_fresnel_schlick/function.h
@@ -20,6 +20,7 @@ class schlick : public nonlinear_function//fresnel
 		schlick()
 		{
 			setParametrization(params::CARTESIAN);
+			setDimX(6);
 		}
 
 		//! \brief Load function specific files
diff --git a/sources/plugins/nonlinear_function_diffuse/function.h b/sources/plugins/nonlinear_function_diffuse/function.h
index ead5217..8025191 100644
--- a/sources/plugins/nonlinear_function_diffuse/function.h
+++ b/sources/plugins/nonlinear_function_diffuse/function.h
@@ -23,6 +23,7 @@ class diffuse_function : public nonlinear_function
         diffuse_function()
         {
             setParametrization(params::CARTESIAN);
+				setDimX(6);
         }
 
 		// Overload the function operator
diff --git a/sources/softs/tests/main.cpp b/sources/softs/tests/main.cpp
index 570e24d..e768d0c 100644
--- a/sources/softs/tests/main.cpp
+++ b/sources/softs/tests/main.cpp
@@ -38,11 +38,61 @@ int main(int argc, char** argv)
 	return nb_tests_failed;
 }
 
+bool is_close(const vec& a, const vec& b)
+{
+	double dist = 0.0;
+	for(int i=0; i<a.size(); ++i)
+	{
+		dist += abs(a[i] - b[i]);
+	}
+	return dist < 1.0E-10;
+}
+
 int parametrization_tests()
 {
 	// Params
 	int nb_tests_failed = 0;
 
+	
+	// Test RUSIN_TH_TD
+	vec rhd(2), cart(6), res(6), vh(3);
+
+	// Convert RUSIN_TH_TD (0,0) to CARTESIAN (0,0,1,0,0,1)
+	rhd[0] = 0;	rhd[1] = 0;
+	res[0] = 0; res[1] = 0; res[2] = 1; res[3] = 0; res[4] = 0; res[5] = 1;
+	params::convert(&rhd[0], params::RUSIN_TH_TD, params::CARTESIAN, &cart[0]);
+	if(!is_close(cart, res)) {
+		nb_tests_failed++;
+		std::cout << "<<FAILED>> out = " << cart << ", while attending " << res << std::endl;
+	}
+
+	// Convert RUSIN_TH_TD (0,pi/2) to CARTESIAN (1,0,0,-1,0,0)
+	rhd[0] = 0;	rhd[1] = M_PI*0.5;
+	res[0] = 1; res[1] = 0; res[2] = 0; res[3] = -1; res[4] = 0; res[5] = 0;
+	params::convert(&rhd[0], params::RUSIN_TH_TD, params::CARTESIAN, &cart[0]);
+	if(!is_close(cart, res)) {
+	  	nb_tests_failed++;
+		std::cout << "<<FAILED>> out = " << cart << ", while attending " << res << std::endl;
+	}
+
+	// Convert RUSIN_TH_TD (0, pi/2) to CARTESIAN (1,0,0,-1,0,0) and to RUSIN_VH (0,0,1)
+	rhd[0] = 0;	rhd[1] = M_PI*0.5;
+	res[0] = 1; res[1] = 0; res[2] = 0; res[3] = -1; res[4] = 0; res[5] = 0;
+	params::convert(&rhd[0], params::RUSIN_TH_TD, params::CARTESIAN, &cart[0]);
+	if(!is_close(cart, res)) {
+	  	nb_tests_failed++;
+		std::cout << "<<FAILED>> out = " << cart << ", while attending " << res << std::endl;
+	}
+	else
+	{
+		params::convert(&rhd[0], params::RUSIN_TH_TD, params::RUSIN_VH, &vh[0]);
+		res.resize(3); res[0] = 0; res[1] = 0; res[2] = 1;
+		if(!is_close(vh, res)) {
+			nb_tests_failed++;
+			std::cout << "<<FAILED>> out = " << vh << ", while attending " << res << std::endl;
+		}
+	}
+/*
 
 	// Test the rotation code
 	vec x(3);
@@ -145,6 +195,6 @@ int parametrization_tests()
     params::convert(&spherical[0], params::SPHERICAL_TL_PL_TV_PV, params::RUSIN_TH_TD_PD, &rusi[0]);
 	 std::cout << "From spherical to rusi" << std::endl;
     std::cout << rusi << std::endl << std::endl;
-
+*/
 	return nb_tests_failed;
 }
-- 
GitLab