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

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

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

55

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

		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) ;
}
		
77
bool rational_fitter_eigen::fit_data(const vertical_segment* d, int np, int nq, rational_function* r) 
78
{
79 80 81 82 83 84 85 86 87 88 89 90
    // 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 ;
91 92 93 94 95 96
}

// 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
97
bool rational_fitter_eigen::fit_data(const vertical_segment* d, int np, int nq, int ny, rational_function_1d* r)
98 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
{
	// 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 ;
			}
		}
	}

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

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

}