rational_fitter.cpp 8.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13
#include "rational_fitter.h"

#include <Eigen/SVD>
#include <Array.hh>
#include <QuadProg++.hh>

#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>

14 15
#include <QTime>

16 17 18 19
#ifdef WIN32
#define isnan(X) ((X != X))
#endif

20 21 22 23 24
fitter* provide_fitter()
{
	return new rational_fitter_quadprog();
}

25 26
data* rational_fitter_quadprog::provide_data() const
{
Laurent Belcour's avatar
Laurent Belcour committed
27
	return new vertical_segment() ;
28 29 30 31 32 33 34
}

function* rational_fitter_quadprog::provide_function() const 
{
	return new rational_function() ;
}

35
rational_fitter_quadprog::rational_fitter_quadprog() : QObject(), _boundary(1.0)
36 37
{
}
38
rational_fitter_quadprog::~rational_fitter_quadprog() 
39 40 41
{
}

42
bool rational_fitter_quadprog::fit_data(const data* dat, function* fit, const arguments &args)
43
{
44
	rational_function* r = dynamic_cast<rational_function*>(fit) ;
Laurent Belcour's avatar
Laurent Belcour committed
45
	const vertical_segment* d = dynamic_cast<const vertical_segment*>(dat) ;
46 47 48 49 50 51 52 53 54 55
	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()) ;
56 57
	r->setMin(d->min()) ;
	r->setMax(d->max()) ;
58 59 60 61

	std::cout << "<<INFO>> np in  [" << _min_np << ", " << _max_np 
	          << "] & nq in [" << _min_nq << ", " << _max_nq << "]" << std::endl ;

62
	int temp_np = _min_np, temp_nq = _min_nq ;
63
	while(temp_np <= _max_np || temp_nq <= _max_nq)
64
	{
65 66
		QTime time ;
		time.start() ;
67
		
68
		if(fit_data(d, temp_np, temp_nq, r))
69
		{
70 71 72 73 74 75 76
			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 ;

77 78 79
			return true ;
		}

80
		std::cout << "<<INFO>> fit using np = " << temp_np << " & nq =  " << temp_nq << " failed\r"  ;
81 82
		std::cout.flush() ;

83
      if(temp_np < _max_np)
84 85 86
		{
			++temp_np ;
		}
87
      if(temp_nq < _max_nq)
88 89 90 91 92 93 94
		{
			++temp_nq ;
		}
	}
	return false ;
}

95
void rational_fitter_quadprog::set_parameters(const arguments& args)
96 97 98 99
{
	_max_np = args.get_float("np", 10) ;
	_max_nq = args.get_float("nq", 10) ;
	_min_np = args.get_float("min-np", _max_np) ;
100
    _min_nq = args.get_float("min-nq", _max_nq) ;
101 102 103 104

    _max_np = std::max<int>(_max_np, _min_np);
    _max_nq = std::max<int>(_max_nq, _min_nq);

105
    _boundary = args.get_float("boundary-constraint", 1.0f);
106 107 108
}
		

Laurent Belcour's avatar
Laurent Belcour committed
109
bool rational_fitter_quadprog::fit_data(const vertical_segment* d, int np, int nq, rational_function* r) 
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
{

	// 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
Laurent Belcour's avatar
Laurent Belcour committed
133
bool rational_fitter_quadprog::fit_data(const vertical_segment* dat, int np, int nq, int ny, rational_function* rf) 
134 135
{
	rational_function* r = dynamic_cast<rational_function*>(rf) ;
Laurent Belcour's avatar
Laurent Belcour committed
136
	const vertical_segment* d = dynamic_cast<const vertical_segment*>(dat) ;
137 138 139 140 141 142
	if(r == NULL || d == NULL)
	{
		std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
		return false ;
	}

143 144 145 146
	// Get the maximum value in data to scale the input parameter space
	// so that it reduces the values of the polynomial
	vec dmax = d->max() ;

147 148
	// Matrices of the problem
	QuadProgPP::Matrix<double> G (0.0, np+nq, np+nq) ;
149
	QuadProgPP::Vector<double> g (0.0, np+nq) ;
150 151
	QuadProgPP::Matrix<double> CI(0.0, np+nq, 2*d->size()) ;
	QuadProgPP::Vector<double> ci(0.0, 2*d->size()) ;
152 153
	QuadProgPP::Matrix<double> CE(0.0, np+nq, 0) ;
	QuadProgPP::Vector<double> ce(0.0, 0) ;
154

155 156
	Eigen::MatrixXd eCI(np+nq, 2*d->size()) ;

157 158 159 160
	// Select the size of the result vector to
	// be equal to the dimension of p + q
	for(int i=0; i<np+nq; ++i)
	{
161
		G[i][i] = 1.0 ; 
162 163 164 165 166 167 168 169 170 171 172
	}
	
	// Each constraint (fitting interval or point
	// add another dimension to the constraint
	// matrix
	for(int i=0; i<d->size(); ++i)	
	{		
		// Norm of the row vector
		double a0_norm = 0.0 ;
		double a1_norm = 0.0 ;

173 174 175
        vec xi = d->get(i) ;

        bool is_boundary = false;
176
        for(int k=0; k<d->dimX(); ++k)
177
        {
178
            is_boundary = is_boundary || (xi[k] <= d->min()[k]) || (xi[k] >= d->max()[k]);
179
        }
Laurent Belcour's avatar
Laurent Belcour committed
180

181 182 183 184 185 186 187 188 189
		// 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 
		// the upper constraint
		for(int j=0; j<np+nq; ++j)
		{
			// Filling the p part
			if(j<np)
			{
190
				const double pi = r->p(xi, j) ;
191 192
				a0_norm += pi*pi ;
				a1_norm += pi*pi ;
193 194
				CI[j][i] =  pi ;
				CI[j][i+d->size()] = -pi ;
195 196 197 198

				// Updating Eigen matrix
				eCI(j,i) = pi ;
				eCI(j,i+d->size()) = -pi ;
199 200 201 202 203 204 205
			}
			// Filling the q part
			else
			{
				vec yl, yu ; 
				d->get(i, yl, yu) ;

206 207 208 209 210 211 212 213
                // Add a constraints for boundary conditions
                if(is_boundary)
                {
                    vec mean = 0.5*(yl+yu);
                    yl = mean + _boundary * (yl - mean);
                    yu = mean + _boundary * (yu - mean);
                }

214
				const double qi = r->q(xi, j-np) ;
215 216
				a0_norm += qi*qi * (yl[ny]*yl[ny]) ;
				CI[j][i] = -yl[ny] * qi ;
217
				
218 219
				a1_norm += qi*qi * (yu[ny]*yu[ny]) ;
				CI[j][i+d->size()] = yu[ny] * qi ;
220 221
				
				// Updating Eigen matrix
222 223
				eCI(j,i) = -yl[ny] * qi ;
				eCI(j,i+d->size()) = yu[ny] * qi ;
224 225 226 227 228
			}
		}
	
		// Set the c vector, will later be updated using the
		// delta parameter.
229 230
		ci[i] = -sqrt(a0_norm) ;
		ci[i+d->size()] = -sqrt(a1_norm) ;
231 232
	}
#ifdef DEBUG
233
	std::cout << "CI = [" ;
234 235 236
	for(int j=0; j<d->size()*2; ++j)
	{
		for(int i=0; i<np+nq; ++i)
237 238 239 240 241 242 243 244
		{
			std::cout << CI[i][j] ;
			if(i != np+nq-1) std::cout << ", ";
		}
		if(j != d->size()*2-1)
			std::cout << ";" << std::endl; 
		else
			std::cout << "]" << std::endl ;
245 246 247 248
	}
#endif
	// Update the ci column with the delta parameter
	// (See Celis et al. 2007 p.12)
249
	Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::HouseholderQRPreconditioner> svd(eCI);
250 251
	const double sigma_m = svd.singularValues()(std::min(2*d->size(), np+nq)-1) ;
	const double sigma_M = svd.singularValues()(0) ;
252

253 254 255 256 257 258 259 260
#ifdef DEBUG
	std::cout << "<<DEBUG>> SVD = [ " ;
	for(int i=0; i<std::min(2*d->size(), np+nq); ++i)
	{
		std::cout << svd.singularValues()(i) << ", " ;
	}
	std::cout << " ]" << std::endl ;
#endif
261
	
262
	double delta = sigma_m / sigma_M ;
263 264

	if(isnan(delta) || (std::abs(delta) == std::numeric_limits<double>::infinity()))
265 266 267 268 269 270 271 272 273 274 275
	{
#ifdef DEBUG
		std::cerr << "<<ERROR>> delta factor is NaN of Inf" << std::endl ;
#endif
		return false ;
	}
	else if(delta == 0.0)
	{
		delta = 1.0 ;
	}

276

Laurent Belcour's avatar
Laurent Belcour committed
277
#ifdef DEBUG
278 279 280 281
	std::cout << "<<DEBUG>> delta factor: " << sigma_m << " / " << sigma_M << " = " << delta << std::endl ;
#endif
	for(int i=0; i<2*d->size(); ++i)	
	{		
282
		ci[i] = ci[i] * delta ; 
283 284 285
#ifdef DEBUG
		std::cout << ci[i] << "\t" ;
#endif
286
	}
287 288
#ifdef DEBUG
	std::cout << std::endl << std::endl ;
289

290 291
	std::cout << eCI << std::endl << std::endl ;
#endif
292 293 294
	// Compute the solution
	QuadProgPP::Vector<double> x;
	double cost = QuadProgPP::solve_quadprog(G, g, CE, ce, CI, ci, x);
295 296


297
	bool solves_qp = !(cost == std::numeric_limits<double>::infinity());
298 299
	for(int i=0; i<np+nq; ++i)
	{
300
		const double v = x[i];
301
		solves_qp = solves_qp && !isnan(v) && (v != std::numeric_limits<double>::infinity()) ;
302 303 304 305 306 307
	}

	if(solves_qp)
	{
		// Recopy the vector d
		std::vector<double> p, q;
308
		double norm = 0.0 ;
309 310
		for(int i=0; i<np+nq; ++i)
		{
311
			const double v = x[i];
312
			norm += v*v ;
313 314
			if(i < np)
			{
315
				p.push_back(v) ;
316 317 318
			}
			else
			{
319
				q.push_back(v) ;
320 321 322
			}
		}

323
		r->update(p, q);
Laurent Belcour's avatar
Laurent Belcour committed
324
#ifdef DEBUG
325
		std::cout << "<<INFO>> got solution " << *r << std::endl ;
Laurent Belcour's avatar
Laurent Belcour committed
326
#endif
327
		return norm > 0.0;
328 329
	}
	else
330

331 332 333 334 335
	{
		return false; 
	}
}

336
Q_EXPORT_PLUGIN2(rational_fitter_quadprog, rational_fitter_quadprog)