Commit 7cc7d399 authored by Laurent Belcour's avatar Laurent Belcour
Browse files

Adding rational fitter using a simple least square fit

parent 233ab13c
......@@ -2,7 +2,9 @@ TEMPLATE = subdirs
SUBDIRS = rational_function \
rational_data \
rational_fitter_cgal \
rational_fitter_quadprog
rational_fitter_quadprog \
rational_fitter_eigen
rational_fitter_cgal.depends = rational_function rational_data
rational_fitter_quadprog.depends = rational_function rational_data
rational_fitter_eigen.depends = rational_function rational_data
......@@ -13,6 +13,8 @@
#include <algorithm>
#include <cmath>
#include <QTime>
typedef CGAL::MP_Float ET ;
typedef CGAL::Quadratic_program<ET> Program ;
typedef CGAL::Quadratic_program_solution<ET> Solution ;
......@@ -55,12 +57,22 @@ bool rational_fitter_cgal::fit_data(const data* dat, function* fit)
std::cout << "<<INFO>> np in [" << _min_np << ", " << _max_np
<< "] & nq in [" << _min_nq << ", " << _max_nq << "]" << std::endl ;
int temp_np = _min_np, temp_nq = _min_nq ;
while(temp_np <= _max_np || temp_nq <= _max_nq)
while(temp_np < _max_np || temp_nq < _max_nq)
{
QTime time ;
time.start() ;
if(fit_data(d, temp_np, temp_nq, r))
{
int msec = time.elapsed() ;
int sec = (msec / 1000) % 60 ;
int min = (msec / 60000) % 60 ;
int hour = (msec / 3600000) ;
std::cout << "<<INFO>> got a fit using np = " << temp_np << " & nq = " << temp_nq << " " << std::endl ;
std::cout << "<<INFO>> it took " << hour << "h " << min << "m " << sec << "s" << std::endl ;
return true ;
}
......
#include "rational_fitter.h"
#include <Eigen/Dense>
#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>
#include <QTime>
rational_fitter_eigen::rational_fitter_eigen() : QObject()
{
}
rational_fitter_eigen::~rational_fitter_eigen()
{
}
bool rational_fitter_eigen::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) ;
if(r == NULL || d == NULL)
{
std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
return false ;
}
// I need to set the dimension of the resulting function to be equal
// to the dimension of my fitting problem
r->setDimX(d->dimX()) ;
r->setDimY(d->dimY()) ;
std::cout << "<<INFO>> np =" << _np << "& nq =" << _nq << std::endl ;
QTime time ;
time.start() ;
if(fit_data(d, _np, _nq, r))
{
int msec = time.elapsed() ;
int sec = (msec / 1000) % 60 ;
int min = (msec / 60000) % 60 ;
int hour = (msec / 3600000) ;
std::cout << "<<INFO>> got a fit" << std::endl ;
std::cout << "<<INFO>> it took " << hour << "h " << min << "m " << sec << "s" << std::endl ;
return true ;
}
std::cout << "<<INFO>> fit failed\r" ;
std::cout.flush() ;
return false ;
}
void rational_fitter_eigen::set_parameters(const arguments& args)
{
_np = args.get_float("np", 10) ;
_nq = args.get_float("nq", 10) ;
}
// TODO Finish
bool rational_fitter_eigen::fit_data(const rational_data* d, int np, int nq, rational_function* r)
{
// Multidimensional coefficients
std::vector<double> Pn ; Pn.reserve(d->dimY()*np) ;
std::vector<double> Qn ; Qn.reserve(d->dimY()*nq) ;
for(int j=0; j<d->dimY(); ++j)
{
if(!fit_data(d, np, nq, j, r))
return false ;
for(int i=0; i<np; ++i) { Pn.push_back(r->getP(i)) ; }
for(int i=0; i<nq; ++i) { Qn.push_back(r->getQ(i)) ; }
}
r->update(Pn, Qn) ;
return true ;
}
// dat is the data object, it contains all the points to fit
// 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_eigen::fit_data(const rational_data* d, int np, int nq, int ny, rational_function* r)
{
// Each constraint (fitting interval or point
// add another dimension to the constraint
// matrix
Eigen::MatrixXd CI(np+nq, d->size()) ;
for(int i=0; i<d->size(); ++i)
{
// 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(int j=0; j<np+nq; ++j)
{
// Filling the p part
if(j<np)
{
const double pi = r->p(d->get(i), j) ;
CI(j, i) = pi ;
}
// Filling the q part
else
{
vec yl, yu ;
d->get(i, yl, yu) ;
const double y = 0.5*(yl[ny] + yu[ny]) ;
const double qi = r->q(d->get(i), j-np) ;
CI(j, i) = -y * qi ;
}
}
}
// Perform the Eigen decomposition of CI CI'
Eigen::MatrixXd M(np+nq, np+nq) ;
M = CI * CI.transpose() ;
Eigen::EigenSolver<Eigen::MatrixXd> solver(M) ;
if(solver.info() == Eigen::Success)
{
// Calculate the minimum eigen value and
// its position
double min_val = solver.eigenvalues()(0).real();
int min_ind = 0;
for(int i=0; i<np+nq; ++i)
{
if(min_val > solver.eigenvalues()(i).real())
{
min_val = solver.eigenvalues()(i).real();
min_ind = i ;
}
}
// Recopy the vector d
std::vector<double> p, q;
p.assign(np, 0.0) ; q.assign(nq, 0.0) ;
for(int i=0; i<np+nq; ++i)
{
if(i < np)
{
p[i] = solver.eigenvectors()(i, min_ind).real() ;
}
else
{
q[i-np] = solver.eigenvectors()(i, min_ind).real() ;
}
}
r->update(p, q) ;
std::cout << "<<INFO>> got solution " << *r << std::endl ;
return true;
}
else
{
return false ;
}
}
Q_EXPORT_PLUGIN2(rational_fitter_eigen, rational_fitter_eigen)
#pragma once
// Include STL
#include <vector>
#include <string>
// Interface
#include <QObject>
#include <core/function.h>
#include <core/data.h>
#include <core/fitter.h>
#include <core/args.h>
#include <rational_function.h>
#include <rational_data.h>
class rational_fitter_eigen : public QObject, public fitter
{
Q_OBJECT
Q_INTERFACES(fitter)
public: // methods
rational_fitter_eigen() ;
virtual ~rational_fitter_eigen() ;
// Fitting a data object
virtual bool fit_data(const data* d, function* fit) ;
// 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 void set_parameters(const arguments& args) ;
protected: // data
// min and Max usable np and nq values for the fitting
int _np, _nq ;
} ;
TARGET = rational_fitter_eigen
TEMPLATE = lib
CONFIG *= qt \
plugin \
eigen
DESTDIR = ../../build
INCLUDEPATH += ../rational_function \
../rational_data \
../..
HEADERS = rational_fitter.h
SOURCES = rational_fitter.cpp
LIBS += -L../../build \
-lrational_function \
-lrational_data
QMAKE_CXXFLAGS += -fPIC
......@@ -11,6 +11,8 @@
#include <algorithm>
#include <cmath>
#include <QTime>
#include <rational_function.h>
#include <rational_data.h>
......@@ -21,18 +23,44 @@ rational_fitter_quadprog::~rational_fitter_quadprog()
{
}
bool rational_fitter_quadprog::fit_data(const data* d, function* fit)
bool rational_fitter_quadprog::fit_data(const data* dat, function* fit)
{
std::cout << "<<INFO>> np in [" << _min_np << ", " << _max_np << "] & nq in [" << _min_nq << ", " << _max_nq << "]" << std::endl ;
rational_function* r = dynamic_cast<rational_function*>(fit) ;
const rational_data* d = dynamic_cast<const rational_data*>(dat) ;
if(r == NULL || d == NULL)
{
std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
return false ;
}
// I need to set the dimension of the resulting function to be equal
// to the dimension of my fitting problem
r->setDimX(d->dimX()) ;
r->setDimY(d->dimY()) ;
std::cout << "<<INFO>> np in [" << _min_np << ", " << _max_np
<< "] & nq in [" << _min_nq << ", " << _max_nq << "]" << std::endl ;
int temp_np = _min_np, temp_nq = _min_nq ;
while(temp_np <= _max_np || temp_nq <= _max_nq)
{
if(fit_data(d, temp_np, temp_nq, fit))
QTime time ;
time.start() ;
if(fit_data(d, temp_np, temp_nq, r))
{
int msec = time.elapsed() ;
int sec = (msec / 1000) % 60 ;
int min = (msec / 60000) % 60 ;
int hour = (msec / 3600000) ;
std::cout << "<<INFO>> got a fit using np = " << temp_np << " & nq = " << temp_nq << " " << std::endl ;
std::cout << "<<INFO>> it took " << hour << "h " << min << "m " << sec << "s" << std::endl ;
return true ;
}
std::cout << "<<INFO>> fitt using np = " << temp_np << " & nq = " << temp_nq << " failed\r" ;
std::cout << "<<INFO>> fit using np = " << temp_np << " & nq = " << temp_nq << " failed\r" ;
std::cout.flush() ;
if(temp_np < _max_np)
......
......@@ -59,25 +59,50 @@ vec rational_function::value(const vec& x) const
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(dimX() == 1)
{
deg[0] = i ;
if(i == 0)
return deg ;
}
int temp_i = i-1 ;
int temp_c ;
while(temp_i > 0)
// 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))
{
temp_c = temp_i % dimX() ;
temp_i = temp_i / dimX() ;
deg[temp_c] += 1 ;
Nk += nk ;
nk *= dimX() ;
++k ;
}
// deg[0] += temp_i ;
// Populate the vector from front to back
int j = i-Nk ;
populate(deg, dimX(), k, j) ;
return deg ;
}
......
......@@ -6,12 +6,30 @@ int main(int argc, int argv)
{
std::ofstream f("input.gnuplot") ;
for(int i=0; i<100; ++i)
const int k = 2 ;
if(k == 1)
{
const float x = i / (float)100.0f ;
const float y = 100.0f * exp(-10.0 * x*x) * x*x - 0.01 *x*x*x ;
f << "#DIM 1 1" << std::endl ;
for(int i=0; i<100; ++i)
{
const float x = i / (float)100.0f ;
const float y = 100.0f * exp(-10.0 * x*x) * x*x - 0.01 *x*x*x ;
f << x << "\t" << y << "\t" << 0.1f << std::endl ;
}
}
else if(k == 2)
{
f << "#DIM 2 1" << std::endl ;
for(int i=0; i<100; ++i)
for(int j=0; j<100; ++j)
{
const float x = i / (float)100.0f ;
const float y = j / (float)100.0f ;
const float z = 100.0f * exp(-10.0 * x*x) * y*y - 0.01 *y*x*y ;
f << x << "\t" << y << "\t" << 0.1f << std::endl ;
f << x << "\t" << y << "\t" << z << "\t" << 0.1f << std::endl ;
}
}
return 0 ;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment