Commit 98488744 authored by Laurent Belcour's avatar Laurent Belcour

Fixing the zero problem comming from the evaluation of degrees

parents faf02c8a 0022d2f4
......@@ -80,9 +80,6 @@ bool compare(std::vector<std::vector<int> > a, std::vector<int> b)
void populate(std::vector<int>& vec, int N, int k, int j)
{
// Already used vectors
std::vector<std::vector<int> > already_used ;
vec[0] = k ;
if(j == 0)
return ;
......@@ -90,8 +87,6 @@ void populate(std::vector<int>& vec, int N, int k, int j)
int tj = j ;
while(tj != 0)
{
already_used.push_back(vec) ;
// First non null index
int nn_index = 0; while(vec[nn_index] == 0) { nn_index = (nn_index+1) % N ; }
......@@ -112,22 +107,40 @@ std::vector<int> rational_function::index2degree(int i) const
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))
if(dimX() == 2)
{
Nk += nk ;
nk *= dimX() ;
++k ;
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
{
// 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) ;
// Populate the vector from front to back
int j = i-Nk ;
populate(deg, dimX(), k, j) ;
}
return deg ;
}
double legendre(double x, int i)
......
TEMPLATE = subdirs
SUBDIRS = \
rational_fitter_cgal \
rational_fitter_quadprog \
rational_fitter_quadproge \
rational_fitter_eigen \
rational_fitter_matlab \
rational_function
rational_fitter_leastsquare \
rational_function \
rational_fitter_matlab
......@@ -54,8 +54,14 @@ bool rational_fitter_cgal::fit_data(const data* dat, function* fit)
r->setDimY(d->dimY()) ;
r->setMin(d->min()) ;
r->setMax(d->max()) ;
/*
for(int i=0; i<20; ++i)
{
std::vector<int> deg = r->index2degree(i) ;
std::cout << deg[0] << ", " << deg[1] << std::endl ;
}
throw ;
*/
std::cout << "<<INFO>> np in [" << _min_np << ", " << _max_np
<< "] & nq in [" << _min_nq << ", " << _max_nq << "]" << std::endl ;
......
#include "rational_fitter.h"
#include <Eigen/Dense>
#include <Eigen/SVD>
#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>
#include <QTime>
data* rational_fitter_leastsquare::provide_data() const
{
return new vertical_segment() ;
}
function* rational_fitter_leastsquare::provide_function() const
{
return new rational_function() ;
}
rational_fitter_leastsquare::rational_fitter_leastsquare() : QObject()
{
}
rational_fitter_leastsquare::~rational_fitter_leastsquare()
{
}
bool rational_fitter_leastsquare::fit_data(const data* dat, function* fit)
{
rational_function* r = dynamic_cast<rational_function*>(fit) ;
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 ;
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()) ;
r->setMin(d->min()) ;
r->setMax(d->max()) ;
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_leastsquare::set_parameters(const arguments& args)
{
_np = args.get_float("np", 10) ;
_nq = args.get_float("nq", 10) ;
_max_iter = args.get_float("max-iter", 1) ;
}
bool rational_fitter_leastsquare::fit_data(const vertical_segment* 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_leastsquare::fit_data(const vertical_segment* d, int np, int nq, int ny, rational_function* r)
{
using namespace Eigen;
double scale = 1;
MatrixXd D(d->size(), np+nq);
VectorXd Y(d->size());
for(int i=0; i<d->size(); ++i)
{
vec yl, yu ;
d->get(i, yl, yu);
const double y = 0.5*(yl[ny] + yu[ny]) ;
Y[i] = y;
vec x = d->get(i);
VectorXd::Map(&x[0], 1) /= scale;
// A row of the constraint matrix has this
// form: [p_{0}(x_i), .., p_{np}(x_i), q_{0}(x_i), .., q_{nq}(x_i)]
for(int j=0; j<np+nq; ++j)
{
if(j<np)
{
const double pi = r->p(x, j) ;
D(i, j) = pi ;
}
else
{
const double qi = r->q(x, j-np) ;
D(i,j) = qi ;
}
}
}
VectorXd pq(np+nq);
// initialize with a constant numerator:
pq.head(np).setZero();
pq(0) = 1;
// Alternate between fixed numerator and fixed denominator
// By default we iterate only once (usually successive iteration introduced more errors)
for(int iter=0; iter<_max_iter; ++iter)
{
// Step 1 fit 1/f_i using 1/q only
{
// Theoretically best weighting scheme to approcimate the true LS problem, which is: sum_i (1/q(x_i) - f_i)^2
MatrixXd A = Y.cwiseAbs2().asDiagonal() * D.rightCols(nq);
VectorXd b = Y.asDiagonal() * (D.leftCols(np) * pq.head(np));
VectorXd x = A.colPivHouseholderQr().solve(b);
// Sometimes, this scheme works better:
// MatrixXd A = Y.asDiagonal() * D.rightCols(nq);
// VectorXd b = D.leftCols(np) * pq.head(np);
// VectorXd x = A.colPivHouseholderQr().solve(b);
pq.tail(nq) = x;
// Statistics on the problem:
// VectorXd sv = A.jacobiSvd().singularValues();
// std::cout << "Singular values: " << sv.transpose() << std::endl;
// std::cout << "Rcond: " << sv(0)/sv(nq-1) << std::endl;
// std::cout << "<<INFO>> LS error " << (A*x-b).norm()/b.norm() << std::endl;
}
VectorXd res = ((D.leftCols(np) * pq.head(np)).array() / (D.rightCols(nq) * pq.tail(nq)).array() - Y.array());
std::cout << "Real LS norm (1): "
<< res.norm() / Y.norm()
<< " ; inf " << res.lpNorm<Infinity>() / Y.lpNorm<Infinity>()
<< " ; L1 " << res.lpNorm<1>() / Y.lpNorm<1>() << std::endl;
// Step 2 fit f_i using p/q with q fix
{
// This scheme gives more weights to small values: (not good)
// VectorXd V = D.rightCols(nq) * pq.tail(nq);
// MatrixXd A = Y.asDiagonal().inverse() * (V.asDiagonal().inverse() * D.leftCols(np));
// VectorXd b = VectorXd::Ones(Y.size());
// VectorXd x = A.colPivHouseholderQr().solve(b);
VectorXd Q = D.rightCols(nq) * pq.tail(nq);
MatrixXd A = Q.asDiagonal().inverse() * D.leftCols(np);
VectorXd x = A.colPivHouseholderQr().solve(Y);
pq.head(np) = x;
// Statistics on the problem:
// VectorXd sv = A.jacobiSvd().singularValues();
// std::cout << "Singular values: " << sv.transpose() << std::endl;
// std::cout << "Rcond: " << sv(0)/sv(np-1) << std::endl;
// std::cout << "<<INFO>> LS error " << (A*x-Y).norm()/Y.norm() << std::endl;
}
} // iterations
std::vector<double> p(np), q(nq);
Eigen::VectorXd::Map(&p[0], np) = pq.head(np);
Eigen::VectorXd::Map(&q[0], nq) = pq.tail(nq);
double s = 1;
for(int i=0; i<np; ++i)
{
p[i] /= s;
s *= scale;
}
s = 1;
for(int i=0; i<nq; ++i)
{
q[i] /= s;
s *= scale;
}
// Evaluate true LS error: sum_i (p(x_i)/q(x_i) - f_i)^2
VectorXd res = ((D.leftCols(np) * pq.head(np)).array() / (D.rightCols(nq) * pq.tail(nq)).array() - Y.array());
std::cout << "<<INFO>> L_2 "
<< res.norm() / Y.norm()
<< " ; L_inf " << res.lpNorm<Infinity>() / Y.lpNorm<Infinity>()
<< " ; L_1 " << res.lpNorm<1>() / Y.lpNorm<1>() << std::endl;
r->update(p, q) ;
std::cout << "<<INFO>> got solution " << *r << std::endl ;
return true;
}
Q_EXPORT_PLUGIN2(rational_fitter_leastsquare, rational_fitter_leastsquare)
#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 <core/rational_function.h>
#include <core/vertical_segment.h>
class rational_fitter_leastsquare : public QObject, public fitter
{
Q_OBJECT
Q_INTERFACES(fitter)
public: // methods
rational_fitter_leastsquare() ;
virtual ~rational_fitter_leastsquare() ;
// Fitting a data object
//
virtual bool fit_data(const data* d, function* fit) ;
// Provide user parameters to the fitter
//
virtual void set_parameters(const arguments& args) ;
// Obtain associated data and functions
//
virtual data* provide_data() const ;
virtual function* provide_function() const ;
protected: // function
// Fitting a data object using np elements in the numerator and nq
// elements in the denominator
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) ;
protected: // data
int _np, _nq ;
int _max_iter;
} ;
TARGET = rational_fitter_leastsquare
TEMPLATE = lib
CONFIG *= qt \
plugin \
eigen
DESTDIR = ../../build
INCLUDEPATH += ../rational_function \
../rational_data \
../..
HEADERS = rational_fitter.h
SOURCES = rational_fitter.cpp
LIBS += -L../../build \
-lcore
#QMAKE_CXXFLAGS += -fPIC
......@@ -214,10 +214,13 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* dat, int np, int
// Set the c vector, will later be updated using the
// delta parameter.
ci(i) = sqrt(a0_norm) ;
ci(i+d->size()) = sqrt(a1_norm) ;
ci(i) = 0.0 ; // sqrt(a0_norm) ;
ci(i+d->size()) = 0.0 ; // sqrt(a1_norm) ;
ce(i) = 0.0 ;
ce(i+d->size()) = 0.0 ;
CI.row(i) /= sqrt(a0_norm) ;
CI.row(i+d->size()) /= sqrt(a1_norm) ;
}
#ifdef DEBUG
std::cout << "CI = [" ;
......@@ -268,7 +271,7 @@ bool rational_fitter_quadprog::fit_data(const vertical_segment* dat, int np, int
#endif
for(int i=0; i<2*d->size(); ++i)
{
CI.row(i) /= ci(i) ;
// CI.row(i) /= ci(i) ;
ci(i) = -delta ;
}
......
......@@ -43,7 +43,8 @@ double rational_function_chebychev::p(const vec& x, int i) const
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
res *= chebychev(x[k], deg[k]);
double xk = 2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5);
res *= chebychev(xk, deg[k]);
}
return res ;
......@@ -54,7 +55,8 @@ double rational_function_chebychev::q(const vec& x, int i) const
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
res *= chebychev(x[k], deg[k]);
double xk = 2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5);
res *= chebychev(xk, deg[k]);
}
return res ;
......
......@@ -128,7 +128,16 @@ int main(int argc, char** argv)
function* f = functions[0] ;
data* d = datas[0] ;
*/
function* f = fit->provide_function() ;
function* f = NULL;
if(args.is_defined("func"))
{
std::cout << "<<INFO>> Using plugin function \"" << args["func"] << "\"" << std::endl ;
f = manager.get_function(args["func"]) ;
}
else
{
f = fit->provide_function() ;
}
data* d = fit->provide_data() ;
d->load(args["input"], args);
......
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