Commit 6a9f5e15 by Laurent Belcour

Not really working Beckmann, need to debug the Jacobian

parent 420c6af9
 ... ... @@ -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 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= 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; ksetMin(_min) ; rs[i]->setMax(_max) ; } return rs[i]; } ... ...
 ... ... @@ -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 rational_function_t : public function virtual void setMin(const vec& min) { function::setMin(min); for(int i=0; isetMin(min); } } virtual void setMax(const vec& max) { function::setMax(max); for(int i=0; isetMax(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; isetMin(min); } } virtual void setMax(const vec& max) { function::setMax(max); for(int i=0; isetMax(max); } } //! \brief Save the rational function to the rational format (see ... ...
 ... ... @@ -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::max()) ; max = args.get_vec("max", _nX, std::numeric_limits::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> v[i] ; vec v = vec::Zero(dimX() + 3*dimY()) ; for(int i=0; i> 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> v[dimX() + i]; v[dimX() + i] /= factor; } for(int i=0; i> v[dimX() + i]; v[dimX() + i] /= factor; } // Check if the data containt a vertical segment around the mean // value. for(int i=0; i> 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 max[i]) // Check if the data containt a vertical segment around the mean // value. for(int i=0; i> 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 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> loaded file \"" << filename << "\"" << std::endl ; ... ...
 ... ... @@ -5,6 +5,7 @@ SConscript('rational_fitter_quadprog/SConscript') # Building functions SConscript('nonlinear_function_diffuse/SConscript') SConscript('nonlinear_function_beckmann/SConscript') SConscript('nonlinear_function_blinn/SConscript') SConscript('nonlinear_function_retroblinn/SConscript') SConscript('rational_function_legendre/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)
 #include "function.h" #include #include #include #include #include #include #include #include ALTA_DLL_EXPORT function* provide_function() { return new beckmann_function(); } // Overload the function operator vec beckmann_function::operator()(const vec& x) const { return value(x); } vec beckmann_function::value(const vec& x) const { vec res(dimY()); double h[3]; params::convert(&x[0], params::CARTESIAN, params::RUSIN_VH, &h[0]); for(int i=0; i 0.0 && x[2]*x[4]>0.0) { res[i] = _ks[i] / (4.0 * x[2]*x[4] * 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; i0.0 && x[2]*x[4]>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[4] * M_PI * a2 * dh2*dh2); // df / dk_s jac[i*nbParameters() + j*2+0] = expo / fac; // df / da_x jac[i*nbParameters() + j*2+1] = - _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> token; if(token.compare("#FUNC") != 0) { std::cerr << "<> parsing the stream. The #FUNC is not the next line defined." << std::endl; return false; } in >> token; if(token.compare("nonlinear_function_beckmann") != 0) { std::cerr << "<> parsing the stream. function name is not the next token." << std::endl; return false; } // Parse the lobe