Commit 11810c0d authored by PACANOWSKI Romain's avatar PACANOWSKI Romain
Browse files

Merge

parents 6bc46350 51c5fa4d
......@@ -56,7 +56,8 @@ NLOPT_DIR = ['#external/build/lib']
NLOPT_LIBS = ['nlopt']
NLOPT_OPT_LIBS = []
# MATLAB library and Engine
## MATLAB library and Engine
##
MATLAB_INC = ['/home/pac/MATLAB/R2013b/extern/include/']
MATLAB_DIR = ['/home/pac/MATLAB/R2013b/bin/glnxa64/']
MATLAB_LIBS = ['']
MATLAB_LIBS = ['engine']
import os
Import('env')
env = env.Clone()
if env.GetOption('clean'):
print "Nothing to do in external"
......
......@@ -127,8 +127,8 @@ void timer::stop()
void timer::reset()
{
_elapsed = 0;
_start = 0;
_stop = 0;
_start = current_time();
_stop = _start;
}
int timer::elapsed() const
......@@ -157,11 +157,8 @@ unsigned int timer::current_time() const
SYSTEMTIME res;
GetSystemTime(&res);
return (unsigned int)(res.wSecond + res.wMinute*60 + res.wHour*360);
#elif __APPLE__
return static_cast<unsigned int>(clock() / CLOCKS_PER_SEC);
#else
struct timespec res;
clock_gettime(CLOCK_REALTIME, &res);
return static_cast<unsigned int>(res.tv_sec);
time_t _t = time(NULL);
return (unsigned int)_t;
#endif
}
......@@ -27,5 +27,4 @@ SOURCES = common.cpp \
vertical_segment.cpp \
rational_function.cpp \
params.cpp \
function.cpp \
# clustering.cpp
function.cpp
......@@ -178,7 +178,25 @@ void params::to_cartesian(const double* invec, params::input intype,
}
break;
case SCHLICK_TL_TK_PROJ_DPHI:
NOT_IMPLEMENTED();
{
// Set the light direction
outvec[3] = sin(invec[0]);
outvec[4] = 0.0;
outvec[5] = cos(invec[0]);
// The view direction is the symmetric of the reflected direction
// with respect to the back direction:
// v = 2 |r.k| k - r
// r = 2 |l.n| n - l = {-lx, -ly, lz }
const double theta = sqrt(invec[1]*invec[1] + invec[2]*invec[2]);
const double Kx = (theta > 0.0) ? (invec[1]/theta)*sin(theta) : 0.0;
const double Ky = (theta > 0.0) ? (invec[2]/theta)*sin(theta) : 0.0;
const double Kz = cos(theta);
const double dotKR = outvec[5]*Kz - outvec[3]*Kx;
outvec[0] = 2.0*dotKR*Kx + outvec[3];
outvec[1] = 2.0*dotKR*Ky;
outvec[2] = 2.0*dotKR*Kz - outvec[5];
}
break;
case RETRO_TL_TVL_PROJ_DPHI:
{
......@@ -187,19 +205,19 @@ void params::to_cartesian(const double* invec, params::input intype,
const double sin_vl = sin(theta_vl);
const double cos_l = cos(invec[0]);
const double sin_l = sin(invec[0]);
const double cos_phi = invec[0] / theta_vl;
const double sin_phi = invec[1] / theta_vl;
const double cos_phi = (theta_vl > 0.0) ? invec[1] / theta_vl : 0.0;
const double sin_phi = (theta_vl > 0.0) ? invec[2] / theta_vl : 0.0;
// Compute the cosine of the outgoing vector using the
// spherical law of cosines
const double cos_v = cos_l*cos_vl + cos_phi*sin_l*sin_vl;
const double sin_v = sqrt(1.0 - cos_v*cos_v);
const double sin_v = sqrt(1.0 - std::min(cos_v*cos_v, 1.0));
if(sin_v != 0.0)
{
const double sin_dphi = (sin_vl*sin_phi) / sin_v;
outvec[0] = sin_v * sqrt(1.0 - sin_dphi*sin_dphi);
outvec[0] = sin_v * sqrt(1.0 - std::min(sin_dphi*sin_dphi, 1.0));
outvec[1] = sin_v * sin_dphi;
}
else
......
......@@ -214,7 +214,8 @@ function* plugins_manager::get_function(const std::string& filename)
// Parameters of the function object
int nX, nY;
params::input param_in; params::output param_out;
params::input param_in = params::UNKNOWN_INPUT;
params::output param_out = params::UNKNOWN_OUTPUT;
arguments args;
// Test for the first line of the file. Should be a ALTA FUNC HEADER
......@@ -331,7 +332,7 @@ function* plugins_manager::get_function(const arguments& args)
}
else
{
std::string filename = args["func"];
const std::string filename = args["func"];
FunctionPrototype myFunction = open_library<FunctionPrototype>(filename, "provide_function");
if(myFunction != NULL)
{
......
......@@ -12,17 +12,14 @@ rational_function_1d::rational_function_1d()
{
}
rational_function_1d::rational_function_1d(int np, int nq, bool separable)
rational_function_1d::rational_function_1d(int nX, unsigned int np, unsigned int nq,
bool separable)
{
resize(np, nq);
_separable = separable;
}
setDimX(nX);
setDimY(1);
rational_function_1d::rational_function_1d(const vec& a,
const vec& b) : a(a), b(b)
{
//update(a, b);
_separable = false;
resize(np, nq);
_separable = separable;
}
bool rational_function_1d::load(std::istream&)
......@@ -34,70 +31,70 @@ void rational_function_1d::save_body(std::ostream& out, const arguments& args) c
{
bool is_matlab = args["export"] == "matlab";
bool is_alta = !args.is_defined("export") || args["export"] == "alta";
if(is_alta)
{
for(int i=0; i<a.size(); ++i)
{
std::vector<int> index = this->index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
out << index[j] << "\t" ;
}
out << a[i] << std::endl ;
}
for(int i=0; i<b.size(); ++i)
{
std::vector<int> index = this->index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
out << index[j] << "\t" ;
}
out << b[i] << std::endl ;
}
}
else if(is_matlab)
{
out << "(";
for(int i=0; i<a.size(); ++i)
{
out << a[i];
std::vector<int> indices = this->index2degree(i);
for(int k=0; k<dimX(); ++k)
{
if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << indices[k];
}
if(i != a.size()-1) { out << " + "; }
}
out << ") / (";
for(int i=0; i<b.size(); ++i)
{
out << b[i] << "x.^" << i;
std::vector<int> indices = this->index2degree(i);
for(int k=0; k<dimX(); ++k)
{
if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << indices[k];
}
if(i != b.size()-1) { out << " + "; }
}
out << ")";
}
else
{
NOT_IMPLEMENTED();
}
const unsigned int np = _p_coeffs.size();
const unsigned int nq = _q_coeffs.size();
if(is_alta)
{
for(unsigned int i=0; i<np; ++i)
{
for(int j=0; j<dimX(); ++j)
{
out << _p_coeffs[i].deg[j] << "\t" ;
}
out << _p_coeffs[i].a << std::endl ;
}
for(unsigned int i=0; i<nq; ++i)
{
for(int j=0; j<dimX(); ++j)
{
out << _q_coeffs[i].deg[j] << "\t" ;
}
out << _q_coeffs[i].a << std::endl ;
}
}
else if(is_matlab)
{
out << "(";
for(unsigned int i=0; i<np; ++i)
{
out << _p_coeffs[i].a;
for(int k=0; k<dimX(); ++k)
{
if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << _q_coeffs[i].deg[k];
}
if(i != np-1) { out << " + "; }
}
out << ") / (";
for(unsigned int i=0; i<nq; ++i)
{
out << _p_coeffs[i].a << "x.^" << i;
for(int k=0; k<dimX(); ++k)
{
if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << _q_coeffs[i].deg[k];
}
if(i != nq-1) { out << " + "; }
}
out << ")";
}
else
{
NOT_IMPLEMENTED();
}
}
void rational_function_1d::update(const vec& in_a,
const vec& in_b)
{
a.resize(in_a.size()) ;
b.resize(in_b.size()) ;
// Get the size of the input vector
const int np = in_a.size();
const int nq = in_b.size();
#define NORMALIZE
......@@ -107,22 +104,45 @@ void rational_function_1d::update(const vec& in_a,
const double b0 = 1.0;
#endif
for(int i=0; i<a.size(); ++i) { a[i] = in_a[i] / b0; }
for(int i=0; i<b.size(); ++i) { b[i] = in_b[i] / b0; }
for(int k=0; k<np; ++k)
{
_p_coeffs[k].a = in_a[k] / b0;
}
for(int k=0; k<nq; ++k)
{
_q_coeffs[k].a = in_b[k] / b0;
}
}
void rational_function_1d::resize(int np, int nq)
void rational_function_1d::resize(unsigned int np, unsigned int nq)
{
const int old_np = a.size();
const int old_nq = b.size();
// It is not possible to resize the multi-dimensional vector
// if the input dimensions size is not defined. This can
// happen at the creation of the rational function object.
if(dimX() == 0) { return; }
// Resize the vector
a.resize(np);
b.resize(nq);
// Resize the numerator
if(_p_coeffs.size() != np)
{
_p_coeffs.resize(np);
for(unsigned int k=0; k<np; ++k)
{
std::vector<int> deg = index2degree(k);
_p_coeffs[k].deg = deg;
}
}
// Set the new coeffs to zero
for(int i=old_np; i<np; ++i) { a[i] = 0.0; }
for(int i=old_nq; i<nq; ++i) { b[i] = 0.0; }
// Resize the denominator
if(_q_coeffs.size() != nq)
{
_q_coeffs.resize(nq);
for(unsigned int k=0; k<nq; ++k)
{
std::vector<int> deg = index2degree(k);
_q_coeffs[k].deg = deg;
}
}
}
......@@ -130,40 +150,26 @@ void rational_function_1d::resize(int np, int nq)
// Get the p_i and q_j function
vec rational_function_1d::p(const vec& x) const
{
vec res(_nY) ;
unsigned int const np = a.size() / _nY ;
vec res(1) ;
for(int k=0; k<_nY; ++k)
const unsigned int np = _p_coeffs.size() ;
for(unsigned int i=0; i<np; ++i)
{
double p = 0.0f ;
for(unsigned int i=0; i<np; ++i)
{
p += a[k*_nY + i]*this->p(x, i) ;
}
res[k] = p ;
res[0] += _p_coeffs[i].a * this->p(x, i) ;
}
return res ;
}
vec rational_function_1d::q(const vec& x) const
{
vec res(_nY) ;
unsigned int const nq = b.size() / _nY ;
vec res(1) ;
for(int k=0; k<_nY; ++k)
const unsigned int np = _q_coeffs.size() ;
for(unsigned int i=0; i<np; ++i)
{
double q = 0.0f ;
for(unsigned int i=0; i<nq; ++i)
{
q += b[k*_nY + i]*this->q(x, i) ;
}
res[k] = q ;
res[0] += _q_coeffs[i].a * this->q(x, i) ;
}
return res ;
}
......@@ -287,19 +293,25 @@ std::vector<int> rational_function_1d::index2degree(int i) const
// Get the p_i and q_j function
double rational_function_1d::p(const vec& x, int i) const
{
std::vector<int> deg = index2degree(i);
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
res *= pow(2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5), deg[k]) ;
//res *= pow(x[k], deg[k]) ;
{
const double xp = 2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5);
res *= pow(xp, _p_coeffs[i].deg[k]) ;
}
return res ;
}
double rational_function_1d::q(const vec& x, int i) const
{
return p(x, i);
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
const double xp = 2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5);
res *= pow(xp, _q_coeffs[i].deg[k]) ;
}
return res ;
}
// Overload the function operator
......@@ -307,20 +319,20 @@ vec rational_function_1d::value(const vec& x) const
{
vec res(1) ;
unsigned int const np = a.size() / _nY ;
unsigned int const nq = b.size() / _nY ;
unsigned int const np = _p_coeffs.size();
unsigned int const nq = _q_coeffs.size();
double p = 0.0f ;
double q = 0.0f ;
for(unsigned int i=0; i<np; ++i)
{
p += a[i]*this->p(x, i) ;
p += _p_coeffs[i].a*this->p(x, i) ;
}
for(unsigned int i=0; i<nq; ++i)
{
q += b[i]*this->q(x, i) ;
q += _q_coeffs[i].a*this->q(x, i) ;
}
res[0] = p/q ;
......@@ -330,25 +342,28 @@ vec rational_function_1d::value(const vec& x) const
std::ostream& operator<< (std::ostream& out, const rational_function_1d& r)
{
const unsigned int np = r._p_coeffs.size();
const unsigned int nq = r._q_coeffs.size();
std::cout << "p = [" ;
for(int i=0; i<r.a.size(); ++i)
for(unsigned int i=0; i<np; ++i)
{
if(i != 0)
{
std::cout << ", " ;
}
std::cout << r.a[i] ;
std::cout << r._p_coeffs[i].a ;
}
std::cout << "]" << std::endl ;
std::cout << "q = [" ;
for(int i=0; i<r.b.size(); ++i)
for(unsigned int i=0; i<nq; ++i)
{
if(i != 0)
{
std::cout << ", " ;
}
std::cout << r.b[i] ;
std::cout << r._q_coeffs[i].a ;
}
std::cout << "]" << std::endl ;
......@@ -384,9 +399,7 @@ rational_function_1d* rational_function::get(int i)
{
if(rs[i] == NULL)
{
rs[i] = new rational_function_1d(np, nq);
rs[i]->setDimX(dimX());
rs[i]->setDimY(dimY());
rs[i] = new rational_function_1d(dimX(), np, nq);
// Test if the input domain is not empty. If one dimension of
// the input domain is a point, I manually inflate this dimension
......
......@@ -21,28 +21,28 @@ class rational_function_1d : public function
public: // methods
rational_function_1d() ;
rational_function_1d(int np, int nq, bool separable = false) ;
rational_function_1d(const vec& a, const vec& b) ;
rational_function_1d(int nX, unsigned int np, unsigned int nq,
bool separable = false) ;
virtual ~rational_function_1d() {}
/* FUNCTION INHERITANCE */
/* FUNCTION INHERITANCE */
// Overload the function operator
//! Overload the function operator
virtual vec value(const vec& x) const ;
virtual vec operator()(const vec& x) const { return value(x) ; }
// IO function to text files
virtual bool load(std::istream& in);
//! IO function to text files
virtual bool load(std::istream& in);
//! \brief Save the rational function expansion. It should
//! not be store in any variable (e.g. "y = rf(x);") as the
//! nD rational function can append factor to the 1D rational
//! function.
virtual void save_body(std::ostream&, const arguments&) const;
//! \brief Save the rational function expansion. It should
//! not be store in any variable (e.g. "y = rf(x);") as the
//! nD rational function can append factor to the 1D rational
//! function.
virtual void save_body(std::ostream&, const arguments&) const;
/* RATIONAL FUNCTION SPECIFIC */
/* RATIONAL FUNCTION SPECIFIC */
//! Evaluate the numerator \f$p(\mathbf{x})\f$ of the rational
//! function. This function is provided to allow fast
......@@ -50,6 +50,7 @@ class rational_function_1d : public function
//! algorithm to efficiently evaluate recursively defined
//! polynomials.
virtual vec p(const vec& x) const ;
//! Evaluate the denominator \f$q(\mathbf{x})\f$ of the rational
//! function. This function is provided to allow fast
//! implementation. For example one can use the Clenshaw
......@@ -64,19 +65,41 @@ class rational_function_1d : public function
//! denominator of the rational function.
virtual double q(const vec& x, int j) const ;
// Update the function
//! Update the coefficient vectors with new values. The new values
//! are normalized by the first element of the denominator
//! coefficients.
virtual void update(const vec& in_a,
const vec& in_b) ;
// Resize the polynomial
virtual void resize(int np, int nq);
//! Resize the polynomial.
virtual void resize(unsigned int np, unsigned int nq);
// Get the coefficients
virtual double getP(int i) const { return a[i]; }
virtual double getQ(int i) const { return b[i]; }
virtual vec getP() const { return a; }
virtual vec getQ() const { return b; }
//! Get the i-th coefficient of the numerator.
virtual double getP(int i) const { return _p_coeffs[i].a; }
//! Get the i-th coefficient of the denominator.
virtual double getQ(int i) const { return _q_coeffs[i].a; }
//! Get the vector of coefficient for the numerator.
virtual vec getP() const
{
const int np = _p_coeffs.size();
vec t(np);
for(int i=0; i<np; ++i) {t[i] = _p_coeffs[i].a; }
return t;
}
//! Get the vector of coefficient for the denominator.
virtual vec getQ() const
{
const int nq = _q_coeffs.size();
vec t(nq);
for(int i=0; i<nq; ++i) {t[i] = _q_coeffs[i].a; }
return t;
}
// STL stream ouput
friend std::ostream& operator<< (std::ostream& out,
......@@ -96,14 +119,28 @@ class rational_function_1d : public function
protected: // data
// Store the coefficients for the moment, I assume
// the functions to be polynomials.
vec a, b ;