Cleaning the Fresnel code to accept multidimensional values.

Corrections in the parsing of arguments
Correcting the Jacobian of Fresnel functions

sources/core/args.h

@@ -30,7 +30,8 @@ class arguments
             for(int i=0; i<argc; ++i)
             {
                 std::string temp(argv[i]) ;
-                std::string key, data ;
+                std::string key;
+                std::string data;
 
                 if(temp.compare(0, 2, "--") == 0)
                 {
@@ -46,10 +47,8 @@ class arguments
                         if(next[0] == '-')
                         {
                             --i;
-                            continue;
                         }
-
-                        if(next[0] == '[' && next[next.size()-1] != ']')
+                        else if(next[0] == '[' && next[next.size()-1] != ']')
                         {
                             data.append(next);
 
@@ -81,8 +80,9 @@ class arguments
 #ifdef DEBUG_ARGS
                     std::cout << "]" << std::endl;
 #endif
+                    _map.insert(std::pair<std::string, std::string>(key, data)) ;
                 }
-                _map.insert(std::pair<std::string, std::string>(key, data)) ;
+
             }
         }
 		~arguments()
@@ -132,7 +132,7 @@ class arguments
 			}
         }
         //! \brief update the value \a val stored under key \a key
-        bool update(const std::string& key, const std::string& val)
+        void update(const std::string& key, const std::string& val)
         {
             _map[key] = val;
         }
@@ -311,6 +311,16 @@ class arguments
 			  return current_args;
 		  }
 
+          friend std::ostream& operator<<(std::ostream& out, const arguments& args)
+          {
+              std::map<std::string, std::string>::const_iterator it;
+              for(it=args._map.begin(); it!=args._map.end(); ++it)
+              {
+                  out<< "[" << it->first << "] -> " << it->second << std::endl;
+              }
+              return out;
+          }
+
 
 	private: // data
 

sources/core/fitter.h

@@ -8,6 +8,7 @@
 /*! \brief Fitting interface for generic fitting algorithms
  *  \ingroup core
  *
+ *  \details
  */
 class fitter
 {

sources/core/function.h

@@ -761,12 +761,12 @@ class fresnel : public nonlinear_function
 			{
 				for(int i=0; i<nb_func_params; ++i)
 				{
-					jac[y + _nY*i] = func_jacobian[y + _nY*i] * fres_value[y];
+                    jac[y*nb_params + _nY*i] = func_jacobian[y*nb_func_params + _nY*i] * fres_value[y];
 				}
 
 				for(int i=0; i<nb_fres_params; ++i)
 				{
-					jac[y + _nY*(i+nb_func_params)] = fres_jacobian[y + _nY*i] * func_value[y];
+                    jac[y*nb_params + _nY*(i+nb_func_params)] = fres_jacobian[y*nb_fres_params + _nY*i] * func_value[y];
 				}
 			}
 

sources/plugins/nonlinear_fitter_ceres/fitter.cpp

@@ -94,6 +94,77 @@ class CeresFunctor : public ceres::CostFunction
 		nonlinear_function* _f;
 };
 
+/*! \brief A CostFunction to fit one color channel at a time.
+ *
+ *  \details This forces the functions to have a full formulation. All the
+ *  parameters need to have a value per color channel.
+ */
+class ColorChannelCost : public ceres::CostFunction
+{
+    public:
+    ColorChannelCost(nonlinear_function* f, const vec xi, int j) : _f(f), _xi(xi), _j(j)
+        {
+            // The function should be set 1D prior to this call
+            assert(f->dimY() == 1);
+
+            set_num_residuals(1);
+            mutable_parameter_block_sizes()->push_back(f->nbParameters());
+        }
+
+        virtual bool Evaluate(double const* const* x, double* y, double** dy) const
+        {
+            // Check that the parameters used are within the bounds defined
+            // by the function
+            vec p_min = _f->getParametersMin();
+            vec p_max = _f->getParametersMax();
+            for(int i=0; i<_f->nbParameters(); ++i)
+            {
+                if(x[0][i] < p_min[i] || x[0][i] > p_max[i])
+                {
+                    return false;
+                }
+            }
+
+            // Update the parameters vector
+            vec _p(_f->nbParameters());
+            for(int i=0; i<_f->nbParameters(); ++i) { _p[i] = x[0][i]; }
+            _f->setParameters(_p);
+
+            double _di = _xi[_f->dimX() + _j];
+
+            // Should add the resulting vector completely
+            y[0] = _di - (*_f)(_xi)[0];
+
+            if(dy != NULL)
+            {
+                df(dy);
+            }
+
+            return true;
+        }
+
+        // The parameter of the function _f should be set prior to this function
+        // call. If not it will produce undesirable results.
+        virtual void df(double ** fjac) const
+        {
+            // Get the jacobian of the function at position x_i for the current
+            // set of parameters (set prior to function call)
+            vec _jac = _f->parametersJacobian(_xi);
+
+            // Fill the columns of the matrix
+            for(int j=0; j<_f->nbParameters(); ++j)
+            {
+                fjac[0][j] = -_jac[j];
+            }
+        }
+
+    protected:
+
+        const vec _xi;
+        int _j;
+        nonlinear_function* _f;
+};
+
 nonlinear_fitter_ceres::nonlinear_fitter_ceres() 
 {
 }
@@ -117,7 +188,7 @@ bool nonlinear_fitter_ceres::fit_data(const data* d, function* fit, const argume
         return false;
     }
     nonlinear_function* nf = dynamic_cast<nonlinear_function*>(fit);
-	 nf->bootstrap(d, args);
+
 
 #ifndef DEBUG
 	 std::cout << "<<DEBUG>> number of parameters: " << nf->nbParameters() << std::endl;
@@ -126,54 +197,113 @@ bool nonlinear_fitter_ceres::fit_data(const data* d, function* fit, const argume
 	 {
 		 return true;
 	 }
+#ifdef FIT_CHANNELS
+     if(args.is_defined("ceres-channels"))
+     {
+         std::cout << "<<WARNING>> will fit the output dimensions separately" << std::endl;
+         std::cout << "<<WARNING>> make sur the function is separable." << std::endl;
+         return fit_channel(d, nf, args);
+     }
+     else
+#endif
+     {
+         /* Bootstrap the function */
+         nf->bootstrap(d, args);
+
+         /* the following starting values provide a rough fit. */
+         vec p = nf->parameters();
+
+         // Create the problem
+         ceres::Problem problem;
+         for(int i=0; i<d->size(); ++i)
+         {
+             vec xi = d->get(i);
+             problem.AddResidualBlock(new CeresFunctor(nf, xi), NULL, &p[0]);
+         }
+
+         // Solves the NL problem
+         ceres::Solver::Summary summary;
+         ceres::Solve(options, &problem, &summary);
+
 
+#ifdef DEBUG
+         std::cout << summary.BriefReport() << std::endl;
+#endif
+         std::cout << "<<INFO>> found parameters: " << p << std::endl;
+
+         nf->setParameters(p);
+         return true;
+     }
+}
+
+#ifdef FIT_CHANNELS
+bool nonlinear_fitter_ceres::fit_channel(const data* d, nonlinear_function* nf,
+                                         const arguments& args)
+{
     /* the following starting values provide a rough fit. */
     vec p = nf->parameters();
 
-	 // Create the problem
-	 ceres::Problem problem;
-	 for(int i=0; i<d->size(); ++i)
-	 {
-		 vec xi = d->get(i);
-		 problem.AddResidualBlock(new CeresFunctor(nf, xi), NULL, &p[0]);
-	 }
+    // Convert the current function to monochromatic
+    nf->setDimY(1);
 
-	 // Solver's options
-	 ceres::Solver::Options options;
-	 if(args.is_defined("ceres-max-num-iterations"))
-	 {
-		 options.max_num_iterations = args.get_int("ceres-max-num-iterations", 50); // Default value = 50
-	 }
-	 
-	 if(args["ceres-factorization"] == "normal-cholesky")
-	 {
-		 options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;
-	 }
-	 else
-	 {
-		 options.linear_solver_type = ceres::DENSE_QR;
-	 }
 
-	 if(args.is_defined("ceres-debug"))
-	 {
-		 options.minimizer_progress_to_stdout = true; // Default value = false;
-	 }
+    for(int c=0; c<d->dimY(); ++c)
+    {
+        nf->bootstrap(d, args);
 
+        // Temp parameter vector
+        vec temp_p = nf->parameters();
 
-     // Solves the NL problem
-     ceres::Solver::Summary summary;
-     ceres::Solve(options, &problem, &summary);
+        // Create the problem
+        ceres::Problem problem;
+        for(int i=0; i<d->size(); ++i)
+        {
+            vec xi = d->get(i);
+            problem.AddResidualBlock(new ColorChannelCost(nf, xi, c), NULL, &temp_p[0]);
+        }
 
+        // Solves the NL problem
+        ceres::Solver::Summary summary;
+        ceres::Solve(options, &problem, &summary);
 
 #ifdef DEBUG
-     std::cout << summary.BriefReport() << std::endl;
+        std::cout << summary.BriefReport() << std::endl;
 #endif
-     std::cout << "<<INFO>> found parameters: " << p << std::endl;
 
-     nf->setParameters(p);
-     return true;
+        // Update the resulting parameter vector
+        for(int i=0; i<nf->nbParameters(); ++i)
+        {
+            p[c*nf->nbParameters() + i] = temp_p[i];
+        }
+    }
+
+    std::cout << "<<INFO>> found parameters: " << p << std::endl;
+    nf->setDimY(d->dimY());
+    nf->setParameters(p);
+
+    return true;
 }
+#endif
+
 
 void nonlinear_fitter_ceres::set_parameters(const arguments& args)
 {
+    if(args.is_defined("ceres-max-num-iterations"))
+    {
+        options.max_num_iterations = args.get_int("ceres-max-num-iterations", 50); // Default value = 50
+    }
+
+    if(args["ceres-factorization"] == "normal-cholesky")
+    {
+        options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;
+    }
+    else
+    {
+        options.linear_solver_type = ceres::DENSE_QR;
+    }
+
+    if(args.is_defined("ceres-debug"))
+    {
+        options.minimizer_progress_to_stdout = true; // Default value = false;
+    }
 }

sources/plugins/nonlinear_fitter_ceres/fitter.h

@@ -11,6 +11,9 @@
 #include <core/args.h>
 #include <core/vertical_segment.h>
 
+// CERES include
+#include <ceres/ceres.h>
+
 /*! \brief A non-linear fitter using the CERES solver
  *  \ingroup plugins
  *
@@ -65,8 +68,16 @@ class nonlinear_fitter_ceres: public fitter
 		virtual void set_parameters(const arguments& args) ;
 
 	protected: // function
-
+#ifdef FIT_CHANNELS
+        // Fit a single color channel from the data
+        //
+        virtual bool fit_channel(const data*, nonlinear_function* fit,
+                                 const arguments& args) ;
+#endif
 
 	protected: // data
+
+        // Fitter options
+        ceres::Solver::Options options;
 } ;
 

sources/plugins/nonlinear_fresnel_schlick/function.cpp

@@ -17,30 +17,33 @@ ALTA_DLL_EXPORT function* provide_function()
 //! Load function specific files
 void schlick::load(std::istream& in)
 {
-	fresnel::load(in);
+    fresnel::load(in);
 
-	// Parse line until the next comment
-	while(in.peek() != '#')
-	{
-		char line[256];
-		in.getline(line, 256);
-	}
+    // Parse line until the next comment
+    while(in.peek() != '#')
+    {
+        char line[256];
+        in.getline(line, 256);
+    }
 
-	// Checking for the comment line #FUNC nonlinear_fresnel_schlick
-	std::string token;
-	in >> token;
-	if(token != "#FUNC") { std::cerr << "<<ERROR>> parsing the stream. The #FUNC is not the next line defined." << std::endl; }
+    // Checking for the comment line #FUNC nonlinear_fresnel_schlick
+    std::string token;
+    in >> token;
+    if(token != "#FUNC") { std::cerr << "<<ERROR>> parsing the stream. The #FUNC is not the next line defined." << std::endl; }
 
-	in >> token;
-	if(token != "nonlinear_fresnel_schlick") { std::cerr << "<<ERROR>> parsing the stream. function name is not the next token." << std::endl; }
+    in >> token;
+    if(token != "nonlinear_fresnel_schlick") { std::cerr << "<<ERROR>> parsing the stream. function name is not the next token." << std::endl; }
 
-	// R [double]
-	in >> token >> R;
+    // R [double]
+    for(int i=0; i<dimY(); ++i)
+    {
+        in >> token >> R[i];
+    }
 }
 
 void schlick::save_call(std::ostream& out, const arguments& args) const
 {
-	bool is_alta   = !args.is_defined("export") || args["export"] == "alta";
+    bool is_alta   = !args.is_defined("export") || args["export"] == "alta";
 
     if(is_alta)
     {
@@ -51,32 +54,41 @@ void schlick::save_call(std::ostream& out, const arguments& args) const
         out << "(";	f->save_call(out, args); out << ")";
     }
 
-	if(is_alta)
-	{
-		out << "#FUNC nonlinear_fresnel_schlick" << std::endl ;
-		out << "R " << R << std::endl;
-		out << std::endl;
-	}
-	else
-	{
-		out << " * schlick_fresnel(L, V, N, X, Y, " << R << ")";
-	}
+    if(is_alta)
+    {
+        out << "#FUNC nonlinear_fresnel_schlick" << std::endl ;
+        for(int i=0; i<dimY(); ++i)
+        {
+            out << "R " << R[i] << std::endl;
+        }
+        out << std::endl;
+    }
+    else
+    {
+        out << " * schlick_fresnel(L, V, N, X, Y, vec3(";
+        for(int i=0; i<dimY(); ++i)
+        {
+            out << R[i];
+            if(i < _nY-1) { out << ", "; }
+        }
+        out << "))";
+    }
 }
 
 void schlick::save_body(std::ostream& out, const arguments& args) const
 {
-	f->save_body(out, args);
-	bool is_shader = args["export"] == "shader" || args["export"] == "explorer";
-
-	if(is_shader)
-	{
-		out << std::endl;
-		out << "vec3 schlick_fresnel(vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y, float R)" << std::endl;
-		out << "{" << std::endl;
-		out << "\tvec3 H = normalize(L + V);" << std::endl;
-		out << "\treturn vec3(R + (1.0f - R) * pow(1.0f - clamp(dot(H,L), 0.0f, 1.0f), 5));" << std::endl;
-		out << "}" << std::endl;
-	}
+    f->save_body(out, args);
+    bool is_shader = args["export"] == "shader" || args["export"] == "explorer";
+
+    if(is_shader)
+    {
+        out << std::endl;
+        out << "vec3 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 R + (vec3(1.0) - R) * pow(1.0f - clamp(dot(H,V), 0.0f, 1.0f), 5);" << std::endl;
+        out << "}" << std::endl;
+    }
 
 }
 
@@ -91,54 +103,74 @@ vec schlick::fresnelValue(const vec& x) const
     vec res(_nY);
     for(int i=0; i<_nY; ++i)
     {
-        res[i] = R + (1.0 - R) * pow(1.0 - clamp(dotVH, 0.0, 1.0), 5.0);
+        res[i] = R[i] + (1.0 - R[i]) * pow(1.0 - clamp(dotVH, 0.0, 1.0), 5.0);
     }
 
     return res;
 }
 
 //! \brief Number of parameters to this non-linear function
-int schlick::nbFresnelParameters() const 
+int schlick::nbFresnelParameters() const
+{
+    return dimY();
+}
+
+vec schlick::getFresnelParametersMin() const
 {
-	return 1;
+    vec m(dimY());
+    for(int i=0; i<dimY(); ++i) { m[i] = 0.0; }
+    return m;
+}
+
+//! Get the vector of min parameters for the function
+vec schlick::getFresnelParametersMax() const
+{
+    vec M(dimY());
+    for(int i=0; i<dimY(); ++i) { M[i] = 1.0; }
+    return M;
 }
 
 //! \brief Get the vector of parameters for the function
-vec schlick::getFresnelParameters() const 
+vec schlick::getFresnelParameters() const
 {
-	vec p(1);
-	p[0] = R;
-	return p;
+    vec p(dimY());
+    for(int i=0; i<dimY(); ++i) { p[i] = R[i]; }
+    return p;
 }
 
 //! \brief Update the vector of parameters for the function
-void schlick::setFresnelParameters(const vec& p) 
+void schlick::setFresnelParameters(const vec& p)
 {
-	R = p[0];
+    for(int i=0; i<dimY(); ++i) { R[i] = p[i]; }
 }
 
 //! \brief Obtain the derivatives of the function with respect to the
-//! parameters. 
-vec schlick::getFresnelParametersJacobian(const vec& x) const 
+//! parameters.
+vec schlick::getFresnelParametersJacobian(const vec& x) const
 {
-	const int nY = dimY();
+    const int nY = dimY();
     double xp[3], cart[6];
     params::convert(&x[0], input_parametrization(), params::RUSIN_VH, xp);
     params::convert(&x[0], input_parametrization(), params::CARTESIAN, cart);
 
     const double dotVH = xp[0]*cart[0] + xp[1]*cart[1] + xp[2]*cart[2];
 
-	vec jac(nY);
-	for(int i=0; i<nY; ++i)
-	{
-        jac[i] = 1.0 - pow(1.0 - clamp(dotVH, 0.0, 1.0), 5.0);
-	}
-
-	return jac;
+    vec jac(nY*nY);
+    for(int i=0; i<nY; ++i)
+        for(int j=0; j<nY; ++j)
+        {
+            if(i == j) {
+                jac[j*dimY() + i] = 1.0 - pow(1.0 - clamp(dotVH, 0.0, 1.0), 5.0);
+            } else {
+                jac[j*dimY() + i] = 0.0;
+            }
+        }
+
+    return jac;
 }
 
 
 void schlick::fresnelBootstrap(const data* d, const arguments& args)
 {
-	R = 0.5;
+    for(int i=0; i<dimY(); ++i) { R[i] = 1.0; }
 }

sources/plugins/nonlinear_fresnel_schlick/function.h

@@ -37,20 +37,10 @@ class schlick : public fresnel
 		virtual void setFresnelParameters(const vec& p) ;
 
 		//! Get the vector of min parameters for the function
-		virtual vec getFresnelParametersMin() const
-		{
-			vec m(1);
-			m[0] = 0.0;
-			return m;
-		}
+        virtual vec getFresnelParametersMin() const;
 
 		//! Get the vector of min parameters for the function
-		virtual vec getFresnelParametersMax() const
-		{
-			vec M(1);
-			M[0] = 1.0;
-			return M;
-		}
+        virtual vec getFresnelParametersMax() const;
 
 		//! \brief Obtain the derivatives of the function with respect to the
 		//! parameters. 
@@ -59,10 +49,16 @@ class schlick : public fresnel
 		//! \brief Boostrap the function by defining the diffuse term
 		virtual void fresnelBootstrap(const data* d, const arguments& args);
 
+        //! \brief resize the parameter vector
+        virtual void setDimY(int nY)
+        {
+            fresnel::setDimY(nY);
+            R.resize(nY);
+        }
 
 	private: // data
 
 		//! Unidimensional Fresnel reflectance at theta = 0
-		double R;
+        vec R;
 } ;
 

sources/plugins/plugins.pro

@@ -15,6 +15,7 @@ SUBDIRS  = \
             nonlinear_fitter_ipopt						\
             nonlinear_fitter_nlopt						\
 				nonlinear_fresnel_schlick					\
+				nonlinear_fresnel_normalized_schlick	\
 				nonlinear_fresnel_retroschlick			\
 				nonlinear_function_diffuse					\
 				nonlinear_function_blinn					\