rational_fitter.cpp 4.75 KB
Newer Older
1 2 3
#include "rational_fitter.h"

#include <Eigen/Dense>
4
#include <Eigen/SVD>
5 6 7 8 9 10 11 12 13 14

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

#include <QTime>

15 16 17
#include <core/common.h>

ALTA_DLL_EXPORT fitter* provide_fitter()
18 19 20 21
{
	return new rational_fitter_eigen();
}

22 23
data* rational_fitter_eigen::provide_data() const
{
24
	return new vertical_segment() ;
25 26 27 28 29 30 31
}

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

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

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

	std::cout << "<<INFO>> np =" << _np << "& nq =" << _nq  << std::endl ;


	QTime time ;
	time.start() ;

62

63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
	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) ;
}
		
87
bool rational_fitter_eigen::fit_data(const vertical_segment* d, int np, int nq, rational_function* r) 
88
{
89 90 91 92 93 94 95 96 97 98 99 100 101 102
    // For each output dimension (color channel for BRDFs) perform
    // a separate fit on the y-1D rational function.
    for(int j=0; j<d->dimY(); ++j)
    {
        rational_function_1d* rs = r->get(j);
        rs->resize(np, nq);

        if(!fit_data(d, np, nq, j, rs))
        {
            return false ;
        }
    }

    return true ;
103 104 105 106 107 108
}

// 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
109
bool rational_fitter_eigen::fit_data(const 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 196 197
		
		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)