rational_fitter.cpp 8.22 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
20
21
#ifdef WIN32
#define isnan(X) ((X != X))
#endif

using namespace std;

22
23
24
25
26
fitter* provide_fitter()
{
	return new rational_fitter_quadprog();
}

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

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

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

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

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

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

79
80
81
			return true ;
		}

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

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

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

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

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

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

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

145
146
147
148
	// 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() ;

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

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

159
160
161
162
	// Select the size of the result vector to
	// be equal to the dimension of p + q
	for(int i=0; i<np+nq; ++i)
	{
163
		G[i][i] = 1.0 ; 
164
165
166
167
168
169
170
171
172
173
174
	}
	
	// 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 ;

175
176
177
        vec xi = d->get(i) ;

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

183
184
185
186
187
188
189
190
191
		// 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)
			{
192
				const double pi = r->p(xi, j) ;
193
194
				a0_norm += pi*pi ;
				a1_norm += pi*pi ;
195
196
				CI[j][i] =  pi ;
				CI[j][i+d->size()] = -pi ;
197
198
199
200

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

208
209
210
211
212
213
214
215
                // 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);
                }

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

255
256
257
258
259
260
261
262
#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
263
	
264
	double delta = sigma_m / sigma_M ;
265
266

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

278

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

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


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

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

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

333
334
335
336
337
	{
		return false; 
	}
}

338
Q_EXPORT_PLUGIN2(rational_fitter_quadprog, rational_fitter_quadprog)