Commit 27538f91 authored by Laurent Belcour's avatar Laurent Belcour

Putting vertical segments and rational functions into the core directory

parent 57b8a973
......@@ -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
#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(x[k], 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(x[k], 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 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 ;
private: // 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 ;
private: // data
// Store the coefficients for the moment, I assume
// the functions to be polynomials.
std::vector<double> a ;
std::vector<double> b ;
} ;
#include "vertical_segment.h"
#include <string>
#include <sstream>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>
void vertical_segment::load(const std::string& filename)
{
arguments args ;
load(filename, args) ;
}
void vertical_segment::load(const std::string& filename, const arguments& args)
{
std::ifstream file(filename.c_str()) ;
if(!file.is_open())
{
std::cerr << "<<ERROR>> unable to open file \"" << filename << "\"" << std::endl ;
}
double min, max ;
min = args.get_float("min", -std::numeric_limits<double>::max()) ;
max = args.get_float("max", std::numeric_limits<double>::max()) ;
_nX = 0 ; _nY = 0 ;
std::vector<int> vs ; int current_vs = 0 ;
double x, y, dy ;
while(file.good())
{
std::string line ;
std::getline(file, line) ;
std::stringstream linestream(line) ;
// Discard incorrect lines
if(linestream.peek() == '#')
{
linestream.ignore(1) ;
std::string comment ;
linestream >> comment ;
if(comment == std::string("DIM"))
{
linestream >> _nX >> _nY ;
vs.assign(dimY(), 0) ;
for(int k=0; k<dimY(); ++k)
{
vs[k] = 0 ;
}
_min.resize(dimX()) ;
_max.resize(dimX()) ;
for(int k=0; k<dimX(); ++k)
{
_min[k] = std::numeric_limits<double>::max() ;
_max[k] = -std::numeric_limits<double>::max() ;
}
}
else if(comment == std::string("VS"))
{
int t ;
linestream >> t ;
vs[current_vs] = t ; ++current_vs ;
}
continue ;
}
else if(line.empty())
{
continue ;
}
vec v ;
v.resize(dimX() + 3*dimY()) ;
for(int i=0; i<dimX(); ++i)
linestream >> v[i] ;
for(int i=0; i<dimY(); ++i)
linestream >> v[dimX() + i] ;
for(int i=0; i<dimY(); ++i)
{
// TODO, the firts case does not account for the
// dimension of the ouput vector
if(vs[i] == 2)
{
linestream >> v[dimX() + dimY()+i] ;
linestream >> v[dimX() + 2*dimY()+i] ;
}
else if(vs[i] == 1)
{
double dt ;
linestream >> dt ;
v[dimX() + dimY()+i] = v[dimX() + i] * (1.0f - dt) ;
v[dimX() + 2*dimY()+i] = v[dimX() + i] * (1.0f + dt) ;
}
else
{
// TODO Specify the delta in case
// Handle multiple dim
float dt = args.get_float("dt", 0.1);
v[dimX() + dimY()+i] = v[dimX() + i] * (1.0f - dt) ;
v[dimX() + 2*dimY()+i] = v[dimX() + i] * (1.0f + dt) ;
}
}
// If data is not in the interval of fit
// TODO: Update to more dims
bool is_in = true ;
for(int i=0; i<dimX(); ++i)
{
if(v[i] < min || v[i] > max)
{
is_in = false ;
}
}
if(!is_in)
{
continue ;
}
_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]) ;
}
}
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 vertical_segment::get(int i, double& x, double& yl, double& yu) const
{
if(i >= (int)_data.size())
{
return false ;
}
x = _data[i][0] ;
yl = _data[i][1] ;
yu = _data[i][2] ;
return true ;
}
void vertical_segment::get(int i, vec& yl, vec& yu) const
{
yl.resize(dimY()) ; yu.resize(dimY()) ;
for(int j=0; j<dimY(); ++j)
{
yl[j] = _data[i][dimX() + dimY() + j] ;
yu[j] = _data[i][dimX() + 2*dimY() + j] ;
}
}
const vec& vertical_segment::operator[](int i) const
{
return _data[i] ;
}
const vec& vertical_segment::get(int i) const
{
return _data[i] ;
}
int vertical_segment::size() const
{
return _data.size() ;
}
vec vertical_segment::min() const
{
return _min ;
}
vec vertical_segment::max() const
{
return _max ;
}
#pragma once
// Include STL
#include <vector>
#include <string>
// Interface
#include "function.h"
#include "data.h"
#include "fitter.h"
#include "args.h"
class vertical_segment : public data
{
public: // methods
// Load data from a file
virtual void load(const std::string& filename) ;
virtual void load(const std::string& filename, const arguments& args) ;
// Acces to data
virtual bool get(int i, double& x, double& yl, double& yu) const ;
virtual const vec& get(int i) const ;
virtual void get(int i, vec& yl, vec& yu) const ;
virtual const vec& operator[](int i) const ;
// Get data size
virtual int size() const ;
// Get min and max input parameters
virtual vec min() const ;
virtual vec max() const ;
private: // data
// Store for each point of data, the upper
// and lower value
std::vector<vec> _data ;
// Store the min and max value on the input
// domain
vec _min, _max ;
} ;
......@@ -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 ;
......@@ -99,7 +99,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 +123,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) ;
......
......@@ -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>
class rational_fitter_cgal : public QObject, public fitter
{
......@@ -41,8 +41,8 @@ class rational_fitter_cgal : public QObject, public fitter
// Fitting a data object using np elements in the numerator and nq
// elements in the denominator
virtual bool fit_data(const rational_data* d, int np, int nq, rational_function* fit) ;
virtual bool fit_data(const rational_data* dat, int np, int nq, int ny, rational_function* fit) ;
virtual bool fit_data(const vertical_segment* d, int np, int nq, rational_function* fit) ;
virtual bool fit_data(const vertical_segment* dat, int np, int nq, int ny, rational_function* fit) ;