Merge

sources/core/rational_function.cpp

@@ -12,16 +12,18 @@ rational_function_1d::rational_function_1d()
 {
 }
 
-rational_function_1d::rational_function_1d(int np, int nq) 
+rational_function_1d::rational_function_1d(int np, int nq, bool separable) 
 {
 	a.resize(np);
 	b.resize(nq);
+	_separable = separable;
 }
 
 rational_function_1d::rational_function_1d(const vec& a, 
                                            const vec& b) :
 	a(a), b(b)
 {
+	_separable = false;
 }
 
 bool rational_function_1d::load(std::istream&)
@@ -162,40 +164,57 @@ std::vector<int> rational_function_1d::index2degree(int i) const
 
 	if(i == 0)
 		return deg ;
-
+	
+	// The case of one dimensional signals is trivial
 	if(dimX() == 1)
 	{
 		deg[0] = i;
+		return deg;
 	}
-	else if(dimX() == 2)
+
+	// Case of two or more dimension signal, we differ the treatment of
+	// separable signals and non separable signals.
+	if(_separable)
 	{
-		int Nk = 1 ;
-		int k  = 1 ;
-		while(!(i >= Nk && i < Nk+k+1))
+		const int index = i-1;
+		const int base  = index / dimX() + 1;
+		for(int k=0; k<dimX(); ++k)
 		{
-			Nk += k+1 ;
-			++k ;
+			deg[k] = (index % dimX() == k) ? base : 0;
 		}
-
-		int r = i-Nk ;
-		deg[0] = k-r;
-		deg[1] = r;
 	}
 	else
 	{
-		int Nk = 1 ;
-		int k  = 1 ;
-		int dk = estimate_dk(k, dimX()) ;
-		while(!(i >= Nk && i < Nk+dk))
+		if(dimX() == 2)
 		{
-			Nk += dk ;
-			++k ;
-			dk = estimate_dk(k, dimX()) ;
+			int Nk = 1 ;
+			int k  = 1 ;
+			while(!(i >= Nk && i < Nk+k+1))
+			{
+				Nk += k+1 ;
+				++k ;
+			}
+
+			int r = i-Nk ;
+			deg[0] = k-r;
+			deg[1] = r;
 		}
+		else
+		{
+			int Nk = 1 ;
+			int k  = 1 ;
+			int dk = estimate_dk(k, dimX()) ;
+			while(!(i >= Nk && i < Nk+dk))
+			{
+				Nk += dk ;
+				++k ;
+				dk = estimate_dk(k, dimX()) ;
+			}
 
-		// Populate the vector from front to back
-		int j = i-Nk ;
-		populate(deg, dimX(), k, j) ;
+			// Populate the vector from front to back
+			int j = i-Nk ;
+			populate(deg, dimX(), k, j) ;
+		}
 	}
 
 	return deg ;
@@ -305,8 +324,23 @@ rational_function_1d* rational_function::get(int i)
 			rs[i] = new rational_function_1d(np, nq);
 			rs[i]->setDimX(dimX());
 			rs[i]->setDimY(dimY());
-			rs[i]->setMin(getMin()) ;
-			rs[i]->setMax(getMax()) ;
+
+			// Test if the input domain is not empty. If one dimension of
+			// the input domain is a point, I manually inflate this dimension
+			// to avoid numerical issues.
+			vec _min = getMin();
+			vec _max = getMax();
+			for(int k=0; k<dimX(); ++k)
+			{
+				if(_min[k] == _max[k]) 
+				{
+					_min[k] -= 0.1;
+					_max[k] += 0.1;
+				}
+			}
+
+			rs[i]->setMin(_min) ;
+			rs[i]->setMax(_max) ;
 		}
 		return rs[i];
 	}

sources/core/rational_function.h

@@ -21,7 +21,7 @@ class rational_function_1d : public function
 	public: // methods
 
 		rational_function_1d() ;
-		rational_function_1d(int np, int nq) ;
+		rational_function_1d(int np, int nq, bool separable = true) ;
         rational_function_1d(const vec& a, const vec& b) ;
 		virtual ~rational_function_1d() {}
 
@@ -75,6 +75,11 @@ class rational_function_1d : public function
 		// Store the coefficients for the moment, I assume
 		// the functions to be polynomials.
         vec a, b ;
+
+		  //! Is the function separable with respect to its input dimensions?
+		  //! \todo Make possible to have only part of the dimensions
+		  //! separable.
+		  bool _separable;
 } ;
 
 #ifdef TODO
@@ -285,19 +290,11 @@ template<class T> class rational_function_t : public function
 		virtual void setMin(const vec& min)
 		{
 			function::setMin(min);
-			for(int i=0; i<dimY(); ++i)
-			{
-				get(i)->setMin(min);
-			}
 		}
 
 		virtual void setMax(const vec& max)
 		{
 			function::setMax(max);
-			for(int i=0; i<dimY(); ++i)
-			{
-				get(i)->setMax(max);
-			}
 		}
 
 		//! \brief Save the rational function to the rational format (see \ref 
@@ -409,19 +406,11 @@ class rational_function : public function
 		virtual void setMin(const vec& min)
 		{
 			function::setMin(min);
-			for(int i=0; i<dimY(); ++i)
-			{
-				get(i)->setMin(min);
-			}
 		}
 
 		virtual void setMax(const vec& max)
 		{
 			function::setMax(max);
-			for(int i=0; i<dimY(); ++i)
-			{
-				get(i)->setMax(max);
-			}
 		}
 
 		//! \brief Save the rational function to the rational format (see 

sources/core/vertical_segment.cpp

@@ -35,7 +35,7 @@ void vertical_segment::load(const std::string& filename, const arguments& args)
 		std::string line ;
 		std::getline(file, line) ;
 		std::stringstream linestream(line) ;
-		
+
 		// Discard incorrect lines
 		if(linestream.peek() == '#')
 		{
@@ -56,7 +56,7 @@ void vertical_segment::load(const std::string& filename, const arguments& args)
 
 				_min.resize(dimX()) ;
 				_max.resize(dimX()) ;
-				
+
 				min = args.get_vec("min", _nX, -std::numeric_limits<float>::max()) ;
 				max = args.get_vec("max", _nX,  std::numeric_limits<float>::max()) ;
 
@@ -72,112 +72,112 @@ void vertical_segment::load(const std::string& filename, const arguments& args)
 				linestream >> t ;
 				vs[current_vs] = t ; ++current_vs ;
 			}
-         else if(comment == std::string("PARAM_IN"))
-         {
+			else if(comment == std::string("PARAM_IN"))
+			{
 				std::string param;
 				linestream >> param;
 				_in_param = params::parse_input(param);
-         }
-         else if(comment == std::string("PARAM_OUT"))
-         {
+			}
+			else if(comment == std::string("PARAM_OUT"))
+			{
 				std::string param;
 				linestream >> param;
 				_out_param = params::parse_output(param);
-         }
+			}
 			continue ;
 		} 
 		else if(line.empty())
 		{
 			continue ;
 		}
-        else
-        {
+		else
+		{
 
-		vec v = vec::Zero(dimX() + 3*dimY()) ;
-		for(int i=0; i<dimX(); ++i)
-			linestream >> v[i] ;
+			vec v = vec::Zero(dimX() + 3*dimY()) ;
+			for(int i=0; i<dimX(); ++i)
+				linestream >> v[i] ;
 
 
-		// Correction of the data by 1/cosine(theta_L)
-		double factor = 1.0;
-		if(args.is_defined("data-correct-cosine"))
-		{
-			double cart[6];
-			params::convert(&v[0], input_parametrization(), params::CARTESIAN, cart);
-			factor = 1.0/cart[5];
-		}
-		// End of correction
+			// Correction of the data by 1/cosine(theta_L)
+			double factor = 1.0;
+			if(args.is_defined("data-correct-cosine"))
+			{
+				double cart[6];
+				params::convert(&v[0], input_parametrization(), params::CARTESIAN, cart);
+				factor = 1.0/cart[5];
+			}
+			// End of correction
 
-		for(int i=0; i<dimY(); ++i)
-		{
-			linestream >> v[dimX() + i];
-			v[dimX() + i] /= factor;
-		}
+			for(int i=0; i<dimY(); ++i)
+			{
+				linestream >> v[dimX() + i];
+				v[dimX() + i] /= factor;
+			}
 
-        // Check if the data containt a vertical segment around the mean
-        // value.
-        for(int i=0; i<dimY(); ++i)
-        {
-            double min_dt = 0.0;
-            double max_dt = 0.0;
-
-
-            if(vs[i] == 2)
-            {
-                linestream >> min_dt ;
-                linestream >> max_dt ;
-                min_dt = min_dt-v[dimX()+i];
-                max_dt = max_dt-v[dimX()+i];
-            }
-            else if(vs[i] == 1)
-            {
-                double dt ;
-                linestream >> dt ;
-                min_dt = -dt;
-                max_dt =  dt;
-            }
-            else
-            {
-                double dt = args.get_float("dt", 0.1f);
-                min_dt = -dt;
-                max_dt =  dt;
-            }
-
-            if(args.is_defined("dt-relative"))
-            {
-                v[dimX() + dimY()+i]   = v[dimX() + i] * (1.0 + min_dt) ;
-                v[dimX() + 2*dimY()+i] = v[dimX() + i] * (1.0 + max_dt) ;
-            }
-            else
-            {
-                v[dimX() + dimY()+i]   = v[dimX() + i] + min_dt ;
-                v[dimX() + 2*dimY()+i] = v[dimX() + i] + max_dt ;
-            }
-        }
-
-		// If data is not in the interval of fit
-		bool is_in = true ;
-		for(int i=0; i<dimX(); ++i)
-		{
-			if(v[i] < min[i] || v[i] > max[i])
+			// Check if the data containt a vertical segment around the mean
+			// value.
+			for(int i=0; i<dimY(); ++i)
+			{
+				double min_dt = 0.0;
+				double max_dt = 0.0;
+
+
+				if(vs[i] == 2)
+				{
+					linestream >> min_dt ;
+					linestream >> max_dt ;
+					min_dt = min_dt-v[dimX()+i];
+					max_dt = max_dt-v[dimX()+i];
+				}
+				else if(vs[i] == 1)
+				{
+					double dt ;
+					linestream >> dt ;
+					min_dt = -dt;
+					max_dt =  dt;
+				}
+				else
+				{
+					double dt = args.get_float("dt", 0.1f);
+					min_dt = -dt;
+					max_dt =  dt;
+				}
+
+				if(args.is_defined("dt-relative"))
+				{
+					v[dimX() + dimY()+i]   = v[dimX() + i] * (1.0 + min_dt) ;
+					v[dimX() + 2*dimY()+i] = v[dimX() + i] * (1.0 + max_dt) ;
+				}
+				else
+				{
+					v[dimX() + dimY()+i]   = v[dimX() + i] + min_dt ;
+					v[dimX() + 2*dimY()+i] = v[dimX() + i] + max_dt ;
+				}
+			}
+
+			// If data is not in the interval of fit
+			bool is_in = true ;
+			for(int i=0; i<dimX(); ++i)
 			{
-				is_in = false ;
+				if(v[i] < min[i] || v[i] > max[i])
+				{
+					is_in = false ;
+				}
+			}
+			if(!is_in)
+			{
+				continue ;
 			}
-		}
-		if(!is_in)
-		{
-			continue ;
-		}
 
-		_data.push_back(v) ;
+			_data.push_back(v) ;
 
-		// Update min and max
-		for(int k=0; k<dimX(); ++k)
-		{
-			_min[k] = std::min(_min[k], v[k]) ;
-			_max[k] = std::max(_max[k], v[k]) ;
+			// Update min and max
+			for(int k=0; k<dimX(); ++k)
+			{
+				_min[k] = std::min(_min[k], v[k]) ;
+				_max[k] = std::max(_max[k], v[k]) ;
+			}
 		}
-        }
 	}
 
 	std::cout << "<<INFO>> loaded file \"" << filename << "\"" << std::endl ;

sources/plugins/SConscript

@@ -5,6 +5,8 @@ SConscript('rational_fitter_quadprog/SConscript')
 
 # Building functions
 SConscript('nonlinear_function_diffuse/SConscript')
+SConscript('nonlinear_function_beckmann/SConscript')
+SConscript('nonlinear_function_retrobeckmann/SConscript')
 SConscript('nonlinear_function_blinn/SConscript')
 SConscript('nonlinear_function_retroblinn/SConscript')
 SConscript('rational_function_legendre/SConscript')

sources/plugins/nonlinear_function_beckmann/SConscript

+env = Environment()
+env.Append(CPPPATH = ['../../../external/build/include', '../../'])
+env.Append(LIBPATH = ['../../../external/build/lib', '../../build'])
+
+sources = ['function.cpp']
+libs = ['-lcore']
+env.SharedLibrary('../../build/nonlinear_function_beckmann', sources, LIBS=libs)

sources/plugins/nonlinear_function_beckmann/function.cpp

+#include "function.h"
+
+#include <string>
+#include <sstream>
+#include <iostream>
+#include <fstream>
+#include <limits>
+#include <algorithm>
+#include <cmath>
+
+#include <core/common.h>
+
+ALTA_DLL_EXPORT function* provide_function()
+{
+    return new beckmann_function();
+}
+		
+vec beckmann_function::G(const vec& x) const
+{
+	vec res(dimY());
+
+	for(int i=0; i<dimY(); ++i)
+	{
+		const double cl = x[2] / (_a[i] * sqrt(1 - x[2]*x[2]));
+		const double cv = x[5] / (_a[i] * sqrt(1 - x[5]*x[5]));
+
+		res[i] = 1.0;
+		if(cl < 1.6)
+		{
+			res[i] *= (3.535*cl + 2.181*cl*cl) / (1 + 2.276*cl + 2.577*cl*cl);
+		}
+		if(cv < 1.6)
+		{
+			res[i] *= (3.535*cv + 2.181*cv*cv) / (1 + 2.276*cv + 2.577*cv*cv);
+		}
+	}
+	return res;
+}
+
+// Overload the function operator
+vec beckmann_function::operator()(const vec& x) const 
+{
+	return value(x);
+}
+vec beckmann_function::value(const vec& x) const 
+{
+	double h[3];
+	params::convert(&x[0], params::CARTESIAN, params::RUSIN_VH, &h[0]);
+
+	// Compute the Shadow term to init res
+	vec res = G(x);
+
+	for(int i=0; i<dimY(); ++i)
+	{
+		const double a    = _a[i];
+		const double a2   = a*a;
+		const double dh2  = h[2]*h[2];
+		const double expo = exp((dh2 - 1.0) / (a2 * dh2));
+
+		if(h[2] > 0.0 && x[2]*x[5]>0.0)
+		{
+			res[i] *= _ks[i] / (4.0 * x[2]*x[5] * M_PI * a2 * dh2*dh2) * expo;
+		}
+		else
+		{
+			res[i] *= 0.0; 
+		}
+	}
+	return res;
+}
+
+// Reset the output dimension
+void beckmann_function::setDimY(int nY)
+{
+    _nY = nY ;
+
+    // Update the length of the vectors
+    _a.resize(_nY) ;
+    _ks.resize(_nY) ;
+}
+
+//! Number of parameters to this non-linear function
+int beckmann_function::nbParameters() const 
+{
+	return 2*dimY();
+}
+
+//! Get the vector of parameters for the function
+vec beckmann_function::parameters() const 
+{
+	vec res(2*dimY());
+	for(int i=0; i<dimY(); ++i)
+	{
+		res[i*2 + 0] = _ks[i];
+		res[i*2 + 1] = _a[i];
+	}
+	return res;
+}
+
+//! \brief get the min values for the parameters
+vec beckmann_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;
+}
+
+//! Update the vector of parameters for the function
+void beckmann_function::setParameters(const vec& p) 
+{
+	for(int i=0; i<dimY(); ++i)
+	{
+		_ks[i] = p[i*2 + 0];
+		_a[i]  = p[i*2 + 1];
+	}
+}
+
+//! Obtain the derivatives of the function with respect to the 
+//! parameters
+//! \todo finish. 
+vec beckmann_function::parametersJacobian(const vec& x) const 
+{
+	double h[3];
+	params::convert(&x[0], params::CARTESIAN, params::RUSIN_VH, h);
+
+	// Get the geometry term
+	vec g = G(x);
+
+    vec jac(dimY()*nbParameters());
+	 for(int i=0; i<dimY(); ++i)
+	 {
+		 for(int j=0; j<dimY(); ++j)
+		 {
+			 if(i == j && h[2]>0.0 && x[2]*x[5]>0.0)
+			 {
+				 const double a    = _a[i];
+				 const double a2   = a*a;
+				 const double dh2  = h[2]*h[2];
+				 const double expo = exp((dh2 - 1.0) / (a2 * dh2));
+				 const double fac  = (4.0 * x[2]*x[5] * M_PI * a2 * dh2*dh2);
+
+				 // df / dk_s
+				 jac[i*nbParameters() + j*2+0] = g[i] * expo / fac;
+
+				 // df / da_x
+				 jac[i*nbParameters() + j*2+1] = - g[i] * _ks[i] * (expo/(4.0*x[2]*x[5])) * ((2* a * h[2])/(M_PI*a2*a2*dh2)) * (1 + (dh2 - 1.0)*h[2]/(a2*dh2*h[2]));
+			 }
+			 else
+			 {
+				 jac[i*nbParameters() + j*2+0] = 0.0;
+				 jac[i*nbParameters() + j*2+1] = 0.0;
+			 }
+		 }
+	 }
+
+    return jac;
+}
+		
+void beckmann_function::bootstrap(const data* d, const arguments& args)
+{
+	for(int i=0; i<dimY(); ++i)
+	{
+		_ks[i] = 1.0;
+		_a[i]  = 1.0;
+	}
+}
+
+//! Load function specific files
+bool beckmann_function::load(std::istream& in)
+{
+	// Parse line until the next comment