rational_fitter.cpp 4.49 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

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

13 14 15
#include <core/common.h>

ALTA_DLL_EXPORT fitter* provide_fitter()
16 17 18 19
{
	return new rational_fitter_eigen();
}

20
rational_fitter_eigen::rational_fitter_eigen() 
21 22 23 24 25 26
{
}
rational_fitter_eigen::~rational_fitter_eigen() 
{
}

27
bool rational_fitter_eigen::fit_data(const ptr<data>& dat, ptr<function>& fit, const arguments&)
28
{
29 30
	ptr<rational_function> r = dynamic_pointer_cast<rational_function>(fit) ;
	if(!r)
31
	{
Laurent Belcour's avatar
Laurent Belcour committed
32 33 34
		std::cerr << "<<ERROR>> not passing the correct function object to the fitter" << std::endl ;
		return false ;
	} 
35 36
	
	ptr<vertical_segment> d = dynamic_pointer_cast<vertical_segment>(dat) ;
37
	if(!d)
Laurent Belcour's avatar
Laurent Belcour committed
38 39
	{
		std::cerr << "<<ERROR>> not passing the correct data object to the fitter" << std::endl ;
40 41 42 43 44 45 46
		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()) ;
47 48
	r->setMin(d->min()) ;
	r->setMax(d->max()) ;
49 50

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

53
    timer time ;
54 55
	time.start() ;

56

57 58
	if(fit_data(d, _np, _nq, r))
	{
59
        time.stop();
60
		std::cout << "<<INFO>> got a fit" << std::endl ;
61
        std::cout << "<<INFO>> it took " << time << std::endl ;
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77

		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) ;
}
		
78
bool rational_fitter_eigen::fit_data(const ptr<vertical_segment>& d, int np, int nq, const ptr<rational_function>& r) 
79
{
80 81 82 83 84 85 86 87 88 89 90 91
    // 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);
        if(!fit_data(d, np, nq, j, rs))
        {
            return false ;
        }
    }

    return true ;
92 93 94 95 96 97
}

// 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
98
bool rational_fitter_eigen::fit_data(const ptr<vertical_segment>& d, int np, int nq, int ny, rational_function_1d* r)
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
{
	// 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 ;
			}
		}
	}

129 130
//	std::cout << CI << std::endl << std::endl ;

131
	// Perform the Eigen decomposition of CI CI'
132
	Eigen::MatrixXd M(np+nq, np+nq) ;
133 134 135 136 137
	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) ;
138 139 140 141
/*
	std::cout << M.size() << std::endl << std::endl ;
	std::cout << M << std::endl << std::endl ;
*/
142 143
	if(solver.info() == Eigen::Success)
	{
144 145 146
/*
		std::cout << solver.eigenvalues() << std::endl ;
*/
147 148
		// Calculate the minimum eigen value and
		// its position
149
#ifdef NOT_WORKING
150
		int    min_id = 0;
151
		double min_val = solver.eigenvalues().minCoeff(&min_id);
152 153 154 155 156 157 158 159 160 161 162 163
#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
164 165 166 167 168
		
		if(min_val == std::numeric_limits<double>::max())
		{
    		return false;
		}
169
#endif
170
		// Recopy the vector d
Laurent Belcour's avatar
Laurent Belcour committed
171
		vec p(np), q(nq);
172 173
		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);
174 175 176 177 178 179 180 181 182 183 184
		
		r->update(p, q) ;
		std::cout << "<<INFO>> got solution " << *r << std::endl ;
		return true;
	}
	else
	{
		return false ;
	}

}