Commit 185878c7 authored by Laurent Belcour's avatar Laurent Belcour

The retro-fresnel is working with the new interface.

I need to debug the back projection to see if it works
parent 629523f8
......@@ -384,293 +384,3 @@ class product_function : public nonlinear_function
nonlinear_function *f1, *f2;
};
/*! \brief A Fresnel interface
* \ingroup core
* \todo Change it to a product_function instead
*/
class fresnel : public nonlinear_function
{
public: // methods
// Overload the function operator
virtual vec operator()(const vec& x) const
{
return value(x);
}
virtual vec value(const vec& x) const
{
vec fres = fresnelValue(x);
// Convert input space to the function
vec xf(f->dimX());
params::convert(&x[0], params::CARTESIAN, f->input_parametrization(), &xf[0]);
vec func = f->value(xf);
return product(fres, func);
}
//! Load function specific files
virtual bool load(std::istream& in)
{
if(f != NULL)
{
return f->load(in);
}
else
{
std::cout << "<<ERROR>> trying to load a Fresnel object with no base class" << std::endl;
return false;
}
}
//! \brief Provide a first rough fit of the function.
virtual void bootstrap(const data* d, const arguments& args)
{
fresnelBootstrap(d, args);
f->bootstrap(d, args);
}
//! Set the dimension of the input space of the function
virtual void setDimX(int nX)
{
function::setDimX(nX);
f->setDimX(nX);
}
//! Set the dimension of the output space of the function
virtual void setDimY(int nY)
{
function::setDimY(nY);
f->setDimY(nY);
}
// Acces to the domain of definition of the function
virtual void setMin(const vec& min)
{
function::setMin(min);
f->setMin(min);
}
virtual void setMax(const vec& max)
{
function::setMax(max);
f->setMax(max);
}
//! Number of parameters to this non-linear function
virtual int nbParameters() const
{
return f->nbParameters() + nbFresnelParameters();
}
//! Get the vector of parameters for the function
virtual vec parameters() const
{
int nb_func_params = f->nbParameters();
int nb_fres_params = nbFresnelParameters();
int nb_params = nb_func_params + nb_fres_params;
vec params(nb_params);
vec func_params = f->parameters();
for(int i=0; i<nb_func_params; ++i)
{
params[i] = func_params[i];
}
vec fres_params = getFresnelParameters();
for(int i=0; i<nb_fres_params; ++i)
{
params[i+nb_func_params] = fres_params[i];
}
return params;
}
//! Get the vector of min parameters for the function
virtual vec getParametersMax() const
{
int nb_func_params = f->nbParameters();
int nb_fres_params = nbFresnelParameters();
int nb_params = nb_func_params + nb_fres_params;
vec params(nb_params);
vec func_params = f->getParametersMax();
for(int i=0; i<nb_func_params; ++i)
{
params[i] = func_params[i];
}
vec fres_params = getFresnelParametersMax();
for(int i=0; i<nb_fres_params; ++i)
{
params[i+nb_func_params] = fres_params[i];
}
return params;
}
//! Get the vector of min parameters for the function
virtual vec getParametersMin() const
{
int nb_func_params = f->nbParameters();
int nb_fres_params = nbFresnelParameters();
int nb_params = nb_func_params + nb_fres_params;
vec params(nb_params);
vec func_params = f->getParametersMin();
for(int i=0; i<nb_func_params; ++i)
{
params[i] = func_params[i];
}
vec fres_params = getFresnelParametersMin();
for(int i=0; i<nb_fres_params; ++i)
{
params[i+nb_func_params] = fres_params[i];
}
return params;
}
//! Update the vector of parameters for the function
virtual void setParameters(const vec& p)
{
int nb_func_params = f->nbParameters();
int nb_fres_params = nbFresnelParameters();
vec func_params(nb_func_params);
for(int i=0; i<nb_func_params; ++i)
{
func_params[i] = p[i];
}
f->setParameters(func_params);
vec fres_params(nb_fres_params);
for(int i=0; i<nb_fres_params; ++i)
{
fres_params[i] = p[i+nb_func_params];
}
setFresnelParameters(fres_params);
}
//! \brief Obtain the derivatives of the function with respect to the
//! parameters.
virtual vec parametersJacobian(const vec& x) const
{
int nb_func_params = f->nbParameters();
int nb_fres_params = nbFresnelParameters();
int nb_params = nb_func_params + nb_fres_params;
// Convert the input value x to the input space of the function
vec xf(f->dimX());
params::convert(&x[0], params::CARTESIAN, f->input_parametrization(), &xf[0]);
vec func_jacobian = f->parametersJacobian(xf);
vec fres_jacobian = getFresnelParametersJacobian(x);
vec func_value = f->value(xf);
vec fres_value = fresnelValue(x);
// F = fresnel; f = function
// d(F * f)(x) /dp = F(x) df(x) /dp + f(x) dF(x) / dp
vec jac(nb_params*_nY);
for(int y=0; y<_nY; ++y)
{
for(int i=0; i<nb_func_params; ++i)
{
jac[y*nb_params + i] = func_jacobian[y*nb_func_params + i] * fres_value[y];
}
for(int i=0; i<nb_fres_params; ++i)
{
jac[y*nb_params + (i+nb_func_params)] = fres_jacobian[y*nb_fres_params + i] * func_value[y];
}
}
return jac;
}
//! \brief set the value for the base function
void setBase(nonlinear_function* fin)
{
f = fin;
}
//! \brief provide the input parametrization of the object.
virtual params::input input_parametrization() const
{
return params::CARTESIAN;
}
//! \brief provide the outout parametrization of the object.
virtual params::output output_parametrization() const
{
return f->output_parametrization();
}
//! \brief can set the input parametrization of a non-parametrized
//! object. Print an error if it is already defined.
virtual void setParametrization(params::input new_param)
{
function::setParametrization(new_param);
f->setParametrization(new_param);
}
//! \brief can set the output parametrization of a non-parametrized
//! function. Throw an exception if it tries to erase a previously
//! defined one.
virtual void setParametrization(params::output new_param)
{
function::setParametrization(new_param);
f->setParametrization(new_param);
}
protected: // methods
//! \brief the interface for the Fresnel code
virtual vec fresnelValue(const vec& x) const = 0;
//! Number of parameters to this non-linear function
virtual int nbFresnelParameters() const = 0;
//! Get the vector of parameters for the function
virtual vec getFresnelParameters() const = 0;
//! Get the vector of min parameters for the function
virtual vec getFresnelParametersMin() const
{
vec m(nbFresnelParameters());
for(int i=0; i<nbFresnelParameters(); ++i)
{
m[i] = -std::numeric_limits<double>::max();
}
return m;
}
//! Get the vector of min parameters for the function
virtual vec getFresnelParametersMax() const
{
vec M(nbFresnelParameters());
for(int i=0; i<nbFresnelParameters(); ++i)
{
M[i] = std::numeric_limits<double>::max();
}
return M;
}
//! Update the vector of parameters for the function
virtual void setFresnelParameters(const vec& p) = 0;
//! \brief Obtain the derivatives of the function with respect to the
//! parameters.
virtual vec getFresnelParametersJacobian(const vec& x) const = 0;
//! \brief Boostrap the function by defining the diffuse term
virtual void fresnelBootstrap(const data* d, const arguments& args) = 0;
protected: //data
//! the base object
nonlinear_function* f;
};
......@@ -242,9 +242,18 @@ void params::from_cartesian(const double* invec, params::input outtype,
const double Kz = invec[2]-invec[5];
const double norm = sqrt(Kx*Kx + Ky*Ky + Kz*Kz);
outvec[0] = Kx / norm;
outvec[1] = Ky / norm;
outvec[2] = Kz / norm;
if(norm > 1.0E-10)
{
outvec[0] = Kx / norm;
outvec[1] = Ky / norm;
outvec[2] = Kz / norm;
}
else
{
outvec[0] = 0.0;
outvec[1] = 0.0;
outvec[2] = 1.0;
}
}
break;
......
#ifndef EXR_IO_H_
#define EXR_IO_H_
#include <stdexcept>
#include <ImfRgbaFile.h>
#include <ImfArray.h>
template<typename FType>
class t_EXR_IO
{
public:
static bool LoadEXR(const char *filename,int& W,int& H, FType *& pix,int nC=3)
{
Imf::RgbaInputFile file(filename);
Imath::Box2i dw = file.dataWindow();
W = dw.max.x - dw.min.x + 1;
H = dw.max.y - dw.min.y + 1;
Imf::Array2D<Imf::Rgba> pixels;
pixels.resizeErase(H, W);
file.setFrameBuffer (&pixels[0][0] - dw.min.x - dw.min.y * W, 1, W);
file.readPixels (dw.min.y, dw.max.y);
pix = new FType[W*H*3] ;
switch(nC)
{
case 3: for(int i=0;i<H;++i)
for(int j=0;j<W;++j)
{
pix[3*(j+i*W)+0] = pixels[H-i-1][j].r ;
pix[3*(j+i*W)+1] = pixels[H-i-1][j].g ;
pix[3*(j+i*W)+2] = pixels[H-i-1][j].b ;
}
break ;
case 1: for(int i=0;i<H;++i)
for(int j=0;j<W;++j)
pix[j+i*W] = 0.3*pixels[H-i-1][j].r + 0.59*pixels[H-i-1][j].g + 0.11*pixels[H-i-1][j].b ;
break ;
default:
throw std::runtime_error("Unexpected case in EXR_IO.") ;
}
return true ;
}
static bool SaveEXR(const char *filename,int W,int H, const FType *pix,int nC=3)
{
Imf::Array2D<Imf::Rgba> pixels(H,W);
/* Convert separated channel representation to per pixel representation */
switch(nC)
{
case 3:
for (int row=0;row<H;row++)
for(int i=0;i<W;i++)
{
Imf::Rgba &p = pixels[H-row-1][i];
p.r = pix[3*(i+row*W)+0] ;
p.g = pix[3*(i+row*W)+1] ;
p.b = pix[3*(i+row*W)+2] ;
p.a = 1.0 ;
}
break ;
case 1:
for (int row=0;row<H;row++)
for(int i=0;i<W;i++)
{
Imf::Rgba &p = pixels[H-row-1][i];
p.r = pix[i+row*W] ;
p.g = pix[i+row*W] ;
p.b = pix[i+row*W] ;
p.a = FType(1.0) ;
}
break ;
default:
throw std::runtime_error("Unexpected case in EXR_IO.") ;
}
Imf::RgbaOutputFile file(filename, W, H, Imf::WRITE_RGBA);
file.setFrameBuffer(&pixels[0][0], 1, W);
file.writePixels(H);
return true ;
}
};
typedef t_EXR_IO<float> EXR_IO ;
#endif
......@@ -188,5 +188,5 @@ vec schlick::parametersJacobian(const vec& x) const
void schlick::bootstrap(const data* d, const arguments& args)
{
for(int i=0; i<dimY(); ++i) { R[i] = 0.1; }
for(int i=0; i<dimY(); ++i) { R[i] = 0.5; }
}
......@@ -14,82 +14,98 @@ ALTA_DLL_EXPORT function* provide_function()
return new retro_schlick();
}
vec retro_schlick::fresnelValue(const vec& x) const
retro_schlick::retro_schlick()
{
double xp[3], cart[6];
params::convert(&x[0], input_parametrization(), params::SCHLICK_VK, xp);
params::convert(&x[0], input_parametrization(), params::CARTESIAN, cart);
setParametrization(params::CARTESIAN);
setDimX(6);
}
vec retro_schlick::value(const vec& x) const
{
double xp[3], cart[6];
params::convert(&x[0], input_parametrization(), params::SCHLICK_VK, xp);
params::convert(&x[0], input_parametrization(), params::CARTESIAN, cart);
const double dotRK = -xp[0]*cart[0] -xp[1]*cart[1] + xp[2]*cart[2];
const double dotRK = -xp[0]*cart[0] -xp[1]*cart[1] + xp[2]*cart[2];
vec res(_nY);
for(int i=0; i<_nY; ++i)
{
res[i] = R[i] + (1.0 - R[i]) * pow(1.0 - clamp(dotRK, 0.0, 1.0), 5.0);
res[i] = R[i] + (1.0 - R[i]) * pow(1.0 - clamp(dotRK, 0.0, 1.0), 5.0);
}
return res;
}
//! \brief Number of parameters to this non-linear function
int retro_schlick::nbFresnelParameters() const
int retro_schlick::nbParameters() const
{
return dimY();
return dimY();
}
//! \brief Get the vector of parameters for the function
vec retro_schlick::getFresnelParameters() const
vec retro_schlick::parameters() const
{
vec p(dimY());
for(int i=0; i<dimY(); ++i) { p[i] = R[i]; }
vec p(dimY());
for(int i=0; i<dimY(); ++i) { p[i] = R[i]; }
return p;
}
//! \brief Update the vector of parameters for the function
void retro_schlick::setFresnelParameters(const vec& p)
void retro_schlick::setParameters(const vec& p)
{
for(int i=0; i<dimY(); ++i) { R[i] = p[i]; }
for(int i=0; i<dimY(); ++i) { R[i] = p[i]; }
}
//! \brief Obtain the derivatives of the function with respect to the
//! parameters.
vec retro_schlick::getFresnelParametersJacobian(const vec& x) const
vec retro_schlick::parametersJacobian(const vec& x) const
{
const int nY = dimY();
double xp[3], cart[6];
params::convert(&x[0], input_parametrization(), params::SCHLICK_VK, xp);
params::convert(&x[0], input_parametrization(), params::CARTESIAN, cart);
const double dotRK = -xp[0]*cart[0] -xp[1]*cart[1] + xp[2]*cart[2];
vec jac(nY);
for(int i=0; i<nY; ++i)
for(int j=0; j<nY; ++j)
{
if(i == j)
{
jac[i*nY + j] = 1.0 - pow(1.0 - clamp(dotRK, 0.0, 1.0), 5.0);
}
else
{
jac[i*nY + j] = 0.0;
}
}
const int nY = dimY();
double xp[3], cart[6];
params::convert(&x[0], input_parametrization(), params::SCHLICK_VK, xp);
params::convert(&x[0], input_parametrization(), params::CARTESIAN, cart);
const double dotRK = -xp[0]*cart[0] -xp[1]*cart[1] + xp[2]*cart[2];
vec jac(nbParameters()*nY);
for(int i=0; i<nY; ++i)
for(int j=0; j<nY; ++j)
{
if(i == j)
{
jac[i*nY + j] = 1.0 - pow(1.0 - clamp(dotRK, 0.0, 1.0), 5.0);
}
else
{
jac[i*nY + j] = 0.0;
}
}
return jac;
}
vec retro_schlick::getParametersMin() const
{
vec m = vec::Zero(dimY());
return m;
}
vec retro_schlick::getParametersMax() const
{
vec M(dimY());
for(int i=0; i<dimY(); ++i) { M[i] = 1.0; }
return M;
}
void retro_schlick::fresnelBootstrap(const data* d, const arguments& args)
void retro_schlick::bootstrap(const data* d, const arguments& args)
{
for(int i=0; i<dimY(); ++i) {R[i] = 0.5; }
for(int i=0; i<dimY(); ++i) {R[i] = 0.5; }
}
//! Load function specific files
bool retro_schlick::load(std::istream& in)
{
fresnel::load(in);
// Parse line until the next comment
while(in.peek() != '#')
{
......@@ -105,40 +121,30 @@ bool retro_schlick::load(std::istream& in)
// Checking for the comment line #FUNC nonlinear_fresnel_retro_schlick
std::string token;
in >> token;
if(token != "#FUNC")
{
std::cerr << "<<ERROR>> parsing the stream. The #FUNC is not the next line defined." << std::endl;
return false;
}
if(token != "#FUNC")
{
std::cerr << "<<ERROR>> parsing the stream. The #FUNC is not the next line defined." << std::endl;
return false;
}
in >> token;
if(token != "nonlinear_fresnel_retro_schlick")
{
std::cerr << "<<ERROR>> parsing the stream. function name is not the next token." << std::endl;
return false;
}
if(token != "nonlinear_fresnel_retro_schlick")
{
std::cerr << "<<ERROR>> parsing the stream. function name is not the next token." << std::endl;
return false;
}
// R [double]
for(int i=0; i<dimY(); ++i) {
in >> token >> R[i];
}
return true;
for(int i=0; i<dimY(); ++i) {
in >> token >> R[i];
}