rational_fitter.cpp 4.97 KB
Newer Older
1
2
/* ALTA --- Analysis of Bidirectional Reflectance Distribution Functions

3
   Copyright (C) 2013, 2013, 2014, 2016 Inria
4
5
6
7
8
9
10

   This file is part of ALTA.

   This Source Code Form is subject to the terms of the Mozilla Public
   License, v. 2.0.  If a copy of the MPL was not distributed with this
   file, You can obtain one at http://mozilla.org/MPL/2.0/.  */

11
12
13
#include "rational_fitter.h"

#include <Eigen/Dense>
14
#include <Eigen/SVD>
15
16
17
18
19
20
21
22

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

23
24
#include <core/common.h>

25
26
using namespace alta;

27
ALTA_DLL_EXPORT fitter* provide_fitter()
28
29
30
31
{
	return new rational_fitter_eigen();
}

32
rational_fitter_eigen::rational_fitter_eigen() 
33
34
35
36
37
38
{
}
rational_fitter_eigen::~rational_fitter_eigen() 
{
}

39
bool rational_fitter_eigen::fit_data(const ptr<data>& dat, ptr<function>& fit, const arguments&)
40
{
41
42
	ptr<rational_function> r = dynamic_pointer_cast<rational_function>(fit) ;
	if(!r)
43
	{
Laurent Belcour's avatar
Laurent Belcour committed
44
45
46
		std::cerr << "<<ERROR>> not passing the correct function object to the fitter" << std::endl ;
		return false ;
	} 
47
48
	
	ptr<vertical_segment> d = dynamic_pointer_cast<vertical_segment>(dat) ;
49
50
	if(!d
     || d->confidence_interval_kind() != vertical_segment::ASYMMETRICAL_CONFIDENCE_INTERVAL)
Laurent Belcour's avatar
Laurent Belcour committed
51
52
	{
		std::cerr << "<<ERROR>> not passing the correct data object to the fitter" << std::endl ;
53
54
55
		return false ;
	}

56
57
  // XXX: FIT and D may have different values of dimX() and dimY(), but
  // this is fine: we convert values as needed in operator().
58
59
	r->setMin(d->min());
	r->setMax(d->max());
60
61

	std::cout << "<<INFO>> np =" << _np << "& nq =" << _nq  << std::endl ;
62
	r->setSize(_np, _nq);
63

64
    timer time ;
65
66
	time.start() ;

67

68
69
	if(fit_data(d, _np, _nq, r))
	{
70
        time.stop();
71
		std::cout << "<<INFO>> got a fit" << std::endl ;
72
        std::cout << "<<INFO>> it took " << time << std::endl ;
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

		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) ;
}
		
89
bool rational_fitter_eigen::fit_data(const ptr<vertical_segment>& d, int np, int nq, const ptr<rational_function>& r) 
90
{
91
92
    // For each output dimension (color channel for BRDFs) perform
    // a separate fit on the y-1D rational function.
93
    for(int j=0; j<d->parametrization().dimY(); ++j)
94
95
96
97
98
99
100
101
102
    {
        rational_function_1d* rs = r->get(j);
        if(!fit_data(d, np, nq, j, rs))
        {
            return false ;
        }
    }

    return true ;
103
104
105
}

// dat is the data object, it contains all the points to fit
Romain Pacanowski's avatar
Romain Pacanowski committed
106
// np and nq are the degree of the Rational Polynomial Function to fit to the data
107
// y is the dimension to fit on the y-data (e.g. R, G or B for RGB signals)
Romain Pacanowski's avatar
Romain Pacanowski committed
108
// the function returns a rational BRDF function and a boolean
109
bool rational_fitter_eigen::fit_data(const ptr<vertical_segment>& d, int np, int nq, int ny, rational_function_1d* r)
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
{
	// 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 ;
			}
		}
	}

140
141
//	std::cout << CI << std::endl << std::endl ;

142
	// Perform the Eigen decomposition of CI CI'
143
	Eigen::MatrixXd M(np+nq, np+nq) ;
144
145
146
147
148
	M.noalias() = CI * CI.transpose();
  // Faster alternative for large np+nq (only the lower triangular half gets computed):
  // M.setZero();
  // M.selfadjointView<Lower>().rankUpdate(CI.transpose());
	Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> solver(M) ;
149
150
151
152
/*
	std::cout << M.size() << std::endl << std::endl ;
	std::cout << M << std::endl << std::endl ;
*/
153
154
	if(solver.info() == Eigen::Success)
	{
155
156
157
/*
		std::cout << solver.eigenvalues() << std::endl ;
*/
158
159
		// Calculate the minimum eigen value and
		// its position
160
#ifdef NOT_WORKING
161
		int    min_id = 0;
162
		double min_val = solver.eigenvalues().minCoeff(&min_id);
163
164
165
166
167
168
169
170
171
172
173
174
#else
		int    min_id = 0 ;
		double min_val = std::numeric_limits<double>::max() ;
		for(int i=0; i<solver.eigenvalues().size(); ++i)
		{
			double value = solver.eigenvalues()[i] ;
			if(value >= 0 && value < min_val)
			{
				min_id  = i ;
				min_val = value ;
			}
		}
Laurent Belcour's avatar
Laurent Belcour committed
175
176
177
178
179
		
		if(min_val == std::numeric_limits<double>::max())
		{
    		return false;
		}
180
#endif
181
		// Recopy the vector d
Laurent Belcour's avatar
Laurent Belcour committed
182
		vec p(np), q(nq);
183
184
		Eigen::VectorXd::Map(&p[0], np) = solver.eigenvectors().col(min_id).head(np);
		Eigen::VectorXd::Map(&q[0], nq) = solver.eigenvectors().col(min_id).tail(nq);
185
186
187
188
189
190
191
192
193
194
195
		
		r->update(p, q) ;
		std::cout << "<<INFO>> got solution " << *r << std::endl ;
		return true;
	}
	else
	{
		return false ;
	}

}