function.cpp 3.31 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11
#include "function.h"

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

#define one_pi 0.31830988618

12 13 14 15 16
ALTA_DLL_EXPORT function* provide_function()
{
    return new shifted_gamma_function();
}

17 18 19 20 21 22 23 24 25 26
// Overload the function operator
vec shifted_gamma_function::operator()(const vec& x) const 
{
	return value(x);
}
vec shifted_gamma_function::value(const vec& x) const 
{
	// shading
	vec lv(3); lv[0] = x[0]; lv[1] = x[1]; lv[2] = x[2];
	vec n(3); n[0] = 0.0; n[1] = 0.0; n[2] = 1.0;
27
	vec ev(3); ev[0] = x[3]; ev[1] = x[4]; ev[2] = x[5];
28 29 30 31 32 33 34 35 36
	vec halfVector = normalize(lv + ev);
	
	double v_h = dot(ev, halfVector);
	double n_h = dot(n, halfVector);
	double n_l = dot(n, lv); 
	double inLight = 1.0;
	if (n_l < 0.0) inLight = 0.0;
	double n_v = dot(n, ev);

37
	return one_pi * inLight * (n_l * rho_d + rho_s.cwiseProduct(D(alpha, p, n_h, K_ap)).cwiseProduct(G1(n_l)).cwiseProduct(G1 (n_v)).cwiseProduct(Fresnel(F_0, F_1, v_h)));
38 39 40 41 42
}

//! Number of parameters to this non-linear function
int shifted_gamma_function::nbParameters() const 
{
43
	return 11*dimY();
44 45 46 47 48
}

//! Get the vector of parameters for the function
vec shifted_gamma_function::parameters() const 
{
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
	const int n = dimY();

	vec res(nbParameters());
	for(int i=0; i<n; ++i) {
		res[i +  0*n] = sh_c[i];
		res[i +  1*n] = sh_theta0[i];
		res[i +  2*n] = sh_k[i];
		res[i +  3*n] = sh_lambda[i];
		res[i +  4*n] = p[i];
		res[i +  5*n] = F_0[i];
		res[i +  6*n] = F_1[i];
		res[i +  7*n] = K_ap[i];
		res[i +  8*n] = rho_d[i];
		res[i +  9*n] = rho_s[i];
		res[i + 10*n] = alpha[i];
	}

	return res;
67 68 69
}

//! Update the vector of parameters for the function
70
void shifted_gamma_function::setParameters(const vec& pi) 
71
{
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
	const int n = dimY();
	for(int i=0; i<n; ++i) {
		sh_c[i]      = pi[i +  0*n];
		sh_theta0[i] = pi[i +  1*n];
		sh_k[i]      = pi[i +  2*n];
		sh_lambda[i] = pi[i +  3*n];
		p[i]         = pi[i +  4*n];
		F_0[i]       = pi[i +  5*n];
		F_1[i]       = pi[i +  6*n];
		K_ap[i]      = pi[i +  7*n];
		rho_d[i]     = pi[i +  8*n];
		rho_s[i]     = pi[i +  9*n];
		alpha[i]     = pi[i + 10*n];
	}

87 88 89
}

//! Obtain the derivatives of the function with respect to the 
90 91 92 93 94 95 96 97 98 99 100 101 102
//! parameters. \TODO
vec shifted_gamma_function::parametersJacobian(const vec& x) const {
	
	const int n = dimY();
    vec jac(n*nbParameters());

	 for(int i=0; i<n; ++i) {
		 for(int j=0; j<nbParameters(); ++j) {
		 	jac[i + j*n] = 0.0;
		 }
	 }

    return jac;
103 104 105 106 107
}

		
vec shifted_gamma_function::Fresnel(const vec& F0, const vec& F1, double V_H) const
{
108 109 110 111 112
	vec F(dimY());
	for(int i=0; i<dimY(); ++i) {
		F[i] = F0[i] - V_H*F1[i] + (1.0 - F0[i])*pow(1.0 - V_H, 5.0);
	}
	return F;
113 114 115 116 117 118 119 120
}

vec shifted_gamma_function::D(const vec& _alpha, const vec& _p, 
                              double cos_h, const vec& _K) const
{
	double cos2 = cos_h*cos_h;
	double tan2 = (1.-cos2)/cos2;

121 122 123 124 125 126 127 128 129 130
	vec D(dimY());
	for(int i=0; i<dimY(); ++i) {
		const double ax      = _alpha[i] + tan2/_alpha[i];
		const double exp_pow = exp(-ax) / pow(ax, _p[i]);

		const double P22 = exp_pow;

		D[i] = P22 * (one_pi / (cos2 * cos2)) * _K[i];
	}
	return D;
131 132 133 134
}

vec shifted_gamma_function::G1(double theta) const
{
135 136 137 138 139 140
	vec G1(dimY());
	for(int i=0; i<dimY(); ++i)
	{
		const double exp_shc = exp(sh_c[i] * pow(std::max<double>(acos(theta) - sh_theta0[i],0.), sh_k[i]));
		G1[i] =  1.0 + sh_lambda[i] * (1.0 - exp_shc);
	}
141
	return G1;
142
}