Commit 3115307d authored by pacanows's avatar pacanows

Merge

parents 57b8a973 5cae0438
......@@ -4,11 +4,15 @@ CONFIG *= static \
DESTDIR = ../build
HEADERS = args.h \
common.h \
data.h \
fitter.h \
function.h \
plugins_manager.h
HEADERS = args.h \
common.h \
data.h \
fitter.h \
function.h \
plugins_manager.h \
vertical_segment.h \
rational_function.h
SOURCES = plugins_manager.cpp
SOURCES = plugins_manager.cpp \
vertical_segment.cpp \
rational_function.cpp
......@@ -17,20 +17,38 @@ class function
// Overload the function operator
virtual vec operator()(const vec& x) const = 0 ;
virtual vec value(const vec& x) const = 0 ;
// IO function to text files
virtual void load(const std::string& filename) = 0 ;
virtual void save(const std::string& filename, const arguments& args) const = 0 ;
virtual int dimX() const { return _nX ; }
virtual int dimY() const { return _nY ; }
virtual int dimX() const { return _nX ; }
virtual int dimY() const { return _nY ; }
virtual void setDimX(int nX) { _nX = nX ; }
virtual void setDimY(int nY) { _nY = nY ; }
// Acces to the domain of definition of the function
virtual void setMin(const vec& min)
{
assert(min.size() == _nX) ;
_min = min ;
}
virtual void setMax(const vec& max)
{
assert(max.size() == _nX) ;
_max = max ;
}
virtual vec getMin() const { return _min ; }
virtual vec getMax() const { return _max ; }
virtual void setDimX(int nX) { _nX = nX ; }
virtual void setDimY(int nY) { _nY = nY ; }
protected: //data
protected:
// Dimension of the function
// Dimension of the function & domain of
// definition.
int _nX, _nY ;
vec _min, _max ;
} ;
//Q_DECLARE_INTERFACE(function, "Fitter.Function")
Q_DECLARE_INTERFACE(function, "Fitter.Function")
#include "rational_function.h"
#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>
rational_function::rational_function() : a(), b()
{
}
rational_function::rational_function(const std::vector<double>& a,
const std::vector<double>& b) :
a(a), b(b)
{
}
rational_function::~rational_function()
{
}
void rational_function::update(const std::vector<double>& in_a,
const std::vector<double>& in_b)
{
a.reserve(in_a.size()) ;
b.reserve(in_b.size()) ;
a = in_a ;
b = in_b ;
}
// Overload the function operator
vec rational_function::value(const vec& x) const
{
vec res(_nY) ;
unsigned int const np = a.size() / _nY ;
unsigned int const nq = b.size() / _nY ;
for(int k=0; k<_nY; ++k)
{
double p = 0.0f ;
double q = 0.0f ;
for(unsigned int i=0; i<np; ++i)
{
p += a[k*_nY + i]*this->p(x, i) ;
}
for(unsigned int i=0; i<nq; ++i)
{
q += b[k*_nY + i]*this->q(x, i) ;
}
res[k] = p/q ;
}
return res ;
}
void populate(std::vector<int>& vec, int N, int k, int j)
{
vec[0] = k ;
if(j == 0)
return ;
int tj = j ;
while(tj != 0)
{
// First non null index
int nn_index = 0; while(vec[nn_index] == 0) { nn_index = (nn_index+1) % N ; }
// Index of the place where to append
int ap_index = (nn_index + 1) % N ; while(vec[ap_index] == k) { ap_index = (ap_index+1) % N ; }
vec[nn_index] -= 1;
vec[ap_index] += 1;
--tj;
}
}
std::vector<int> rational_function::index2degree(int i) const
{
std::vector<int> deg ; deg.assign(dimX(), 0) ;
if(i == 0)
return deg ;
// Calculate the power (number of elements to put in
// the vector) at which the index is definite.
int Nk = 1 ;
int nk = dimX() ;
int k = 1 ;
while(!(i >= Nk && i < Nk+nk))
{
Nk += nk ;
nk *= dimX() ;
++k ;
}
// Populate the vector from front to back
int j = i-Nk ;
populate(deg, dimX(), k, j) ;
return deg ;
}
double legendre(double x, int i)
{
if(i == 0)
{
return 1;
}
else if(i == 1)
{
return x;
}
else
{
return ((2*i-1)*x*legendre(x, i-1) - (i-1)*legendre(x, i-2)) / (double)i ;
}
}
//#define POLYNOMIALS
// Get the p_i and q_j function
double rational_function::p(const vec& x, int i) const
{
std::vector<int> deg = index2degree(i);
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
#ifdef POLYNOMIALS
res *= pow(x[k], deg[k]) ;
#else // LEGENDRE
res *= legendre(2.0*((x[k] - _min[k]) / _max[k] - 0.5), deg[k]);
#endif
}
return res ;
}
double rational_function::q(const vec& x, int i) const
{
std::vector<int> deg = index2degree(i);
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
#ifdef POLYNOMIALS
res *= pow(x[k], deg[k]) ;
#else // LEGENDRE
res *= legendre(2.0*((x[k] - _min[k]) / _max[k] - 0.5), deg[k]);
#endif
}
return res ;
}
// IO function to text files
void rational_function::load(const std::string& filename)
{
}
void rational_function::save(const std::string& filename, const arguments& args) const
{
save_rational_function(filename) ;
}
void rational_function::save_gnuplot(const std::string& filename, const data* d, const arguments& args) const
{
std::ofstream file(filename.c_str(), std::ios_base::trunc);
for(int i=0; i<d->size(); ++i)
{
vec v = d->get(i) ;
// vec y1 ; y1.assign(d->dimY(), 0.0) ;
// for(int k=0; k<d->dimY(); ++k) { y1[k] = v[d->dimX() + k] ; }
vec y2 = value(v) ;
for(int u=0; u<d->dimX(); ++u)
file << v[u] << "\t" ;
for(int u=0; u<d->dimY(); ++u)
file << y2[u] << "\t" ;
file << std::endl ;
}
}
void rational_function::save_rational_function(const std::string& filename) const
{
std::ofstream file(filename.c_str(), std::ios_base::trunc);
file << "#DIM " << _nX << " " << _nY << std::endl ;
file << "#NP " << a.size() / _nY << std::endl ;
file << "#NQ " << b.size() / _nY << std::endl ;
file << "#BASIS poly" << std::endl ;
for(unsigned int i=0; i<a.size() / _nY; ++i)
{
std::vector<int> index = index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
file << index[j] << "\t" ;
}
file << a[i] << std::endl ;
}
for(unsigned int i=0; i<b.size(); ++i)
{
std::vector<int> index = index2degree(i) ;
for(unsigned int j=0; j<index.size() / _nY; ++j)
{
file << index[j] << "\t" ;
}
file << b[i] << std::endl ;
}
}
std::ostream& operator<< (std::ostream& out, const rational_function& r)
{
std::cout << "p = [" ;
for(unsigned int i=0; i<r.a.size(); ++i)
{
if(i != 0)
{
std::cout << ", " ;
}
std::cout << r.a[i] ;
}
std::cout << "]" << std::endl ;
std::cout << "q = [" ;
for(unsigned int i=0; i<r.b.size(); ++i)
{
if(i != 0)
{
std::cout << ", " ;
}
std::cout << r.b[i] ;
}
std::cout << "]" << std::endl ;
return out ;
}
//Q_EXPORT_PLUGIN2(rational_function, rational_function)
#pragma once
// Include STL
#include <vector>
#include <string>
// Interface
#include <QObject>
#include "function.h"
#include "data.h"
#include "fitter.h"
#include "args.h"
#include "common.h"
class rational_function : public QObject, public function
{
Q_OBJECT
Q_INTERFACES(function)
public: // methods
rational_function() ;
rational_function(const std::vector<double>& a, const std::vector<double>& b) ;
virtual ~rational_function() ;
// Overload the function operator
virtual vec value(const vec& x) const ;
virtual vec operator()(const vec& x) const { return value(x) ; }
// Get the p_i and q_j function
virtual double p(const vec& x, int i) const ;
virtual double q(const vec& x, int j) const ;
// IO function to text files
virtual void load(const std::string& filename) ;
virtual void save(const std::string& filename, const arguments& args) const ;
// Update the function
virtual void update(const std::vector<double>& in_a,
const std::vector<double>& in_b) ;
// Get the coefficients
virtual double getP(int i) const { return a[i] ; }
virtual double getQ(int i) const { return b[i] ; }
// STL stream ouput
friend std::ostream& operator<< (std::ostream& out, const rational_function& r) ;
// Save the rational function to a given format
void save_rational_function(const std::string& filename) const ;
void save_gnuplot(const std::string& filename, const data* d, const arguments& args) const ;
protected: // functions
// Convert an index in N to a vector of degree for a
// multinomial coeffcient. The resulting vector v should
// be used as prod_k x[k]^v[k]
std::vector<int> index2degree(int i) const ;
protected: // data
// Store the coefficients for the moment, I assume
// the functions to be polynomials.
std::vector<double> a ;
std::vector<double> b ;
} ;
#include "rational_data.h"
#include "vertical_segment.h"
#include <string>
#include <sstream>
......@@ -8,13 +8,13 @@
#include <algorithm>
#include <cmath>
void rational_data::load(const std::string& filename)
void vertical_segment::load(const std::string& filename)
{
arguments args ;
load(filename, args) ;
}
void rational_data::load(const std::string& filename, const arguments& args)
void vertical_segment::load(const std::string& filename, const arguments& args)
{
std::ifstream file(filename.c_str()) ;
if(!file.is_open())
......@@ -70,7 +70,7 @@ void rational_data::load(const std::string& filename, const arguments& args)
}
continue ;
}
else if(line.empty()/*!boost::regex_match(line,e)*/)
else if(line.empty())
{
continue ;
}
......@@ -134,27 +134,13 @@ void rational_data::load(const std::string& filename, const arguments& args)
}
}
//TODO Test for small data
/* std::vector<std::vector<double> > temp ;
std::ofstream temp_out("temp.gnuplot", std::ios_base::trunc) ;
for(int i=0; i<20; ++i)
{
int k = (i * _data.size()) / 20 ;
temp_out << _data[k][0] << "\t" << _data[k][1] << "\t" << 0.5*(_data[k][2] - _data[k][1]) << std::endl ;
temp.push_back(_data[k]) ;
}
_data = temp ;
*/
// Sort the vector
// std::sort(_data.begin(), _data.end(), [](const std::vector<double>& a, const std::vector<double>& b){return (a[0]<b[0]);});
std::cout << "<<INFO>> loaded file \"" << filename << "\"" << std::endl ;
std::cout << "<<INFO>> data inside " << _min << " ... " << _max << std::endl ;
std::cout << "<<INFO>> loading data file of R^" << dimX() << " -> R^" << dimY() << std::endl ;
std::cout << "<<INFO>> " << _data.size() << " elements to fit" << std::endl ;
}
bool rational_data::get(int i, double& x, double& yl, double& yu) const
bool vertical_segment::get(int i, double& x, double& yl, double& yu) const
{
if(i >= (int)_data.size())
{
......@@ -168,7 +154,7 @@ bool rational_data::get(int i, double& x, double& yl, double& yu) const
return true ;
}
void rational_data::get(int i, vec& yl, vec& yu) const
void vertical_segment::get(int i, vec& yl, vec& yu) const
{
yl.resize(dimY()) ; yu.resize(dimY()) ;
for(int j=0; j<dimY(); ++j)
......@@ -178,28 +164,26 @@ void rational_data::get(int i, vec& yl, vec& yu) const
}
}
const vec& rational_data::operator[](int i) const
const vec& vertical_segment::operator[](int i) const
{
return _data[i] ;
}
const vec& rational_data::get(int i) const
const vec& vertical_segment::get(int i) const
{
return _data[i] ;
}
int rational_data::size() const
int vertical_segment::size() const
{
return _data.size() ;
}
vec rational_data::min() const
vec vertical_segment::min() const
{
return _min ;
}
vec rational_data::max() const
vec vertical_segment::max() const
{
return _max ;
}
//Q_EXPORT_PLUGIN2(rational_data, rational_data)
......@@ -5,18 +5,13 @@
#include <string>
// Interface
//#include <QObject>
#include <core/function.h>
#include <core/data.h>
#include <core/fitter.h>
#include <core/args.h>
#include "function.h"
#include "data.h"
#include "fitter.h"
#include "args.h"
class rational_data : /*public QObject,*/ public data
class vertical_segment : public data
{
/*
Q_OBJECT
Q_INTERFACES(data)
*/
public: // methods
// Load data from a file
......
TEMPLATE = subdirs
SUBDIRS = rational_function \
rational_data \
SUBDIRS = \
rational_fitter_cgal \
rational_fitter_quadprog \
rational_fitter_quadproge \
rational_fitter_eigen \
rational_fitter_matlab
rational_fitter_matlab \
rational_function
rational_fitter_cgal.depends = rational_function rational_data
rational_fitter_quadprog.depends = rational_function rational_data
rational_fitter_quadproge.depends = rational_function rational_data
rational_fitter_eigen.depends = rational_function rational_data
rational_fitter_matlab.depends = rational_function rational_data
TEMPLATE = lib
CONFIG *= static \
qt \
DESTDIR = ../../build
INCLUDEPATH += ../..
HEADERS = rational_data.h
SOURCES = rational_data.cpp
#QMAKE_CXXFLAGS += -frounding-math -fPIC
......@@ -23,7 +23,7 @@ typedef CGAL::Quadratic_program_options Options ;
data* rational_fitter_cgal::provide_data() const
{
return new rational_data() ;
return new vertical_segment() ;
}
function* rational_fitter_cgal::provide_function() const
......@@ -41,7 +41,7 @@ rational_fitter_cgal::~rational_fitter_cgal()
bool rational_fitter_cgal::fit_data(const data* dat, function* fit)
{
rational_function* r = dynamic_cast<rational_function*>(fit) ;
const rational_data* d = dynamic_cast<const rational_data*>(dat) ;
const vertical_segment* d = dynamic_cast<const vertical_segment*>(dat) ;
if(r == NULL || d == NULL)
{
std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
......@@ -52,6 +52,8 @@ bool rational_fitter_cgal::fit_data(const data* dat, function* fit)
// to the dimension of my fitting problem
r->setDimX(d->dimX()) ;
r->setDimY(d->dimY()) ;
r->setMin(d->min()) ;
r->setMax(d->max()) ;
std::cout << "<<INFO>> np in [" << _min_np << ", " << _max_np
<< "] & nq in [" << _min_nq << ", " << _max_nq << "]" << std::endl ;
......@@ -99,7 +101,7 @@ void rational_fitter_cgal::set_parameters(const arguments& args)
_min_nq = args.get_float("min-nq", _max_nq) ;
}
bool rational_fitter_cgal::fit_data(const rational_data* d, int np, int nq, rational_function* r)
bool rational_fitter_cgal::fit_data(const vertical_segment* d, int np, int nq, rational_function* r)
{
// Multidimensional coefficients
......@@ -123,7 +125,7 @@ bool rational_fitter_cgal::fit_data(const rational_data* d, int np, int nq, rati
// np and nq are the degree of the RP to fit to the data
// y is the dimension to fit on the y-data (e.g. R, G or B for RGB signals)
// the function return a ration BRDF function and a boolean
bool rational_fitter_cgal::fit_data(const rational_data* d, int np, int nq, int ny, rational_function* r)
bool rational_fitter_cgal::fit_data(const vertical_segment* d, int np, int nq, int ny, rational_function* r)
{
// by default, we have a nonnegative QP with Ax - b >= 0
Program qp (CGAL::LARGER, false, 0, false, 0) ;
......@@ -151,11 +153,7 @@ bool rational_fitter_cgal::fit_data(const rational_data* d, int np, int nq, int
double a1_norm = 0.0 ;
vec xi = d->get(i) ;
for(int k=0; k<d->dimX(); ++k)
{
xi[k] /= dmax[k] ;
}
// A row of the constraint matrix has this
// form: [p_{0}(x_i), .., p_{np}(x_i), -f(x_i) q_{0}(x_i), .., -f(x_i) q_{nq}(x_i)]
// For the lower constraint and negated for
......@@ -318,11 +316,11 @@ bool rational_fitter_cgal::fit_data(const rational_data* d, int np, int nq, int
if(i < np)
{
p[i] = v / r->p(dmax, i) ;
p[i] = v ;
}
else
{
q[i-np] = v / r->q(dmax, i-np) ;
q[i-np] = v ;
}
}
r->update(p, q) ;
......
......@@ -7,12 +7,12 @@
// Interface
#include <QObject>
#include <core/function.h>
#include <core/rational_function.h>
#include <core/data.h>
#include <core/vertical_segment.h>
#include <core/fitter.h>
#include <core/args.h>
#include <rational_function.h>
#include <rational_data.h>