Adding nonlinear function retro Beckmann

sources/plugins/SConscript

@@ -6,6 +6,7 @@ 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_retrobeckmann/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_retrobeckmann', sources, LIBS=libs)

sources/plugins/nonlinear_function_retrobeckmann/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::SCHLICK_VK, &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::SCHLICK_VK, 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
+	while(in.peek() != '#')
+	{
+		char line[256];
+		in.getline(line, 256);
+
+		// If we cross the end of the file, or the badbit is
+		// set, the file cannot be loaded
+		if(!in.good())
+			return false;
+	}
+
+    // Checking for the comment line #FUNC nonlinear_function_lafortune
+	std::string token;
+	in >> token;
+	if(token.compare("#FUNC") != 0) 
+	{ 
+		std::cerr << "<<ERROR>> parsing the stream. The #FUNC is not the next line defined." << std::endl; 
+        return false;
+	}
+
+	in >> token;
+   if(token.compare("nonlinear_function_retrobeckmann") != 0) 
+	{
+		std::cerr << "<<ERROR>> parsing the stream. function name is not the next token." << std::endl; 
+        return false;
+	}
+
+	// Parse the lobe
+	for(int i=0; i<_nY; ++i)
+	{
+
+		in >> token >> _ks[i];
+		in >> token >> _a[i];
+	}
+    return true;
+}
+
+
+void beckmann_function::save_call(std::ostream& out, const arguments& args) const
+{
+    bool is_alta   = !args.is_defined("export") || args["export"] == "alta";
+
+    if(is_alta)
+    {
+		out << "#FUNC nonlinear_function_retrobeckmann" << std::endl ;
+
+		 for(int i=0; i<_nY; ++i)
+		 {
+			 out << "Ks " << _ks[i] << std::endl;
+			 out << "a  " << _a[i]  << std::endl;
+		 }
+	
+		 out << std::endl;
+	 }
+	 else
+	 {
+		 out << "retrobeckmann(L, V, N, X, Y, vec3(";
+		 for(int i=0; i<_nY; ++i)
+		 {
+			 out << _ks[i];
+			 if(i < _nY-1) { out << ", "; }
+		 }
+
+		 out << "), vec3(";
+		 for(int i=0; i<_nY; ++i)
+		 {
+			 out << _a[i];
+			 if(i < _nY-1) { out << ", "; }
+		 }
+		 out << "))";
+	 }
+
+}
+
+void beckmann_function::save_body(std::ostream& out, const arguments& args) const
+{
+    bool is_shader = args["export"] == "shader" || args["export"] == "explorer";
+
+    if(is_shader)
+    {
+        out << "vec3 g_beckmann(vec3 M, vec3 N, vec3 a)" << std::endl;
+		  out << "{" << std::endl;
+        out << "\tfloat d = dot(M,N);" << std::endl;
+		  out << "\tvec3 c = d / (a * sqrt(1.0f-d));" << std::endl;
+		  out << "\tvec3 r;" << std::endl;
+		  out << "\tif(c.x < 1.6f) {" << std::endl;
+		  out << "\t\tr.x = (3.535*c.x + 2.181*c.x*c.x) / (1 + 2.276*c.x + 2.577*c.x*c.x);" << std::endl;
+		  out << "\t} else {" << std::endl;
+		  out << "\t\tr.x = 1.0f;" << std::endl;
+		  out << "\t}" << std::endl;
+		  out << "\tif(c.y < 1.6f) {" << std::endl;
+		  out << "\t\tr.y = (3.535*c.y + 2.181*c.y*c.y) / (1 + 2.276*c.y + 2.577*c.y*c.y);" << std::endl;
+		  out << "\t} else {" << std::endl;
+		  out << "\t\tr.y = 1.0f;" << std::endl;
+		  out << "\t}" << std::endl;
+		  out << "\tif(c.z < 1.6f) {" << std::endl;
+		  out << "\t\tr.z = (3.535*c.z + 2.181*c.z*c.z) / (1 + 2.276*c.z + 2.577*c.z*c.z);" << std::endl;
+		  out << "\t} else {" << std::endl;
+		  out << "\t\tr.z = 1.0f;" << std::endl;
+		  out << "\t}" << std::endl;
+		  out << "\treturn r;" << std::endl;
+        out << "}" << std::endl;
+		  out << std::endl;
+        out << "vec3 retrobeckmann(vec3 L, vec3 V, vec3 N, vec3 X, vec3 Y, vec3 ks, vec3 a)" << std::endl;
+        out << "{" << std::endl;
+        out << "\tvec3  R   = 2*dot(V,N)*N - V;" << std::endl;
+        out << "\tvec3  B   = normalize(L + R);" << std::endl;
+        out << "\tfloat bn  = dot(B,N);" << std::endl;
+		  out << "\tfloat ln  = dot(L,N);" << std::endl;
+		  out << "\tfloat vn  = dot(V,N);" << std::endl;
+		  out << "\t" << std::endl;
+        out << "\treturn ks / (4 * " << M_PI << " * a*a * ln*vn) * exp((bn*bn - 1.0) / (a*a*bn*bn)) * g_beckmann(L,N,a) * g_beckmann(V,N,a);" << std::endl;
+        out << "}" << std::endl;
+    }
+}

sources/plugins/nonlinear_function_retrobeckmann/function.h

+#pragma once
+
+// Include STL
+#include <vector>
+#include <string>
+
+// Interface
+#include <core/function.h>
+#include <core/data.h>
+#include <core/fitter.h>
+#include <core/args.h>
+#include <core/common.h>
+
+/*! \brief A Microfacet model using a Gaussian distribution for the
+ *  micro-gemeotry distribution parametrized by the half-angle.
+ *
+ *  \details
+ *
+ *  <h3>Plugin parameters</h3>
+ *
+ */
+class beckmann_function : public nonlinear_function
+{
+
+	public: // methods
+
+		beckmann_function()
+		{
+			setParametrization(params::CARTESIAN);
+			setDimX(6);
+		}
+
+		// Overload the function operator
+		virtual vec operator()(const vec& x) const ;
+		virtual vec value(const vec& x) const ;
+
+		// Geometrical term of the microfacets
+		virtual vec G(const vec& x) const;
+
+		//! \brief Load function specific files
+		virtual bool load(std::istream& in) ;
+
+		//! \brief Export function
+		virtual void save_call(std::ostream& out, const arguments& args) const;
+		virtual void save_body(std::ostream& out, const arguments& args) const;
+
+		//! \brief Boostrap the function by defining the diffuse term
+		//!
+		//! \details
+		virtual void bootstrap(const data* d, const arguments& args);
+
+		//! \brief Number of parameters to this non-linear function
+		virtual int nbParameters() const ;
+
+		//! \brief Get the vector of parameters for the function
+		virtual vec parameters() const ;
+
+		//! \brief get the min values for the parameters
+		virtual vec getParametersMin() const;
+
+		//! \brief Update the vector of parameters for the function
+		virtual void setParameters(const vec& p) ;
+
+		//! \brief Obtain the derivatives of the function with respect to the
+		//! parameters. 
+		virtual vec parametersJacobian(const vec& x) const ;
+
+		//! \brief Provide the dimension of the input space of the function
+		virtual int dimX() const
+		{
+			return 6;
+		}
+
+		//! \brief Set the number of output dimensions
+		void setDimY(int nY);
+
+	private: // data
+
+		vec _ks, _a; // Lobes data
+} ;
+

sources/plugins/nonlinear_function_retrobeckmann/nonlinear_function_retrobeckmann.pro

+TEMPLATE        = lib
+CONFIG         *= plugin \
+                  eigen
+
+DESTDIR         = ../../build
+
+INCLUDEPATH    += ../..
+HEADERS         = function.h
+SOURCES         = function.cpp
+
+
+LIBS           += -L../../build  \
+						-lcore
+

sources/plugins/plugins.pro

@@ -19,6 +19,7 @@ SUBDIRS  = \
 				nonlinear_fresnel_retroschlick			\
 				nonlinear_function_diffuse					\
 				nonlinear_function_beckmann				\
+				nonlinear_function_retrobeckmann			\
 				nonlinear_function_blinn					\
 				nonlinear_function_retroblinn				\
 				nonlinear_function_ward 					\