Commit 46083a70 authored by Laurent Belcour's avatar Laurent Belcour

[Refactoring] Moving the interpolation plugins to the directory

data_interpolants.
parent 56584bd7
# # Building data plugins
#SConscript('data_merl/SConscript')
SConscript('data_interpolant/SConscript')
SConscript('data_interpolants/SConscript')
SConscript('data_io/SConscript')
......
Import('env')
env = env.Clone()
env.AppendUnique(LIBS = env['PLUGIN_LIB'])
build_lib = False
test_env = env
test_env.AppendUnique(LIBS = env['FLANN_LIB'])
test_env.AppendUnique(LIBPATH = env['FLANN_DIR'])
test_env.AppendUnique(CPPPATH = env['FLANN_INC'])
conf = Configure(test_env)
if conf.CheckLibWithHeader(env['FLANN_LIB'], 'flann/flann.hpp', 'CXX'):
build_lib = True
env.AppendUnique(LIBS = env['FLANN_LIB'])
env.AppendUnique(LIBPATH = env['FLANN_DIR'])
env.AppendUnique(CPPPATH = env['FLANN_INC'])
#end
test_env = conf.Finish()
test_env = env
test_env.AppendUnique(LIBS = env['CGAL_LIB'])
test_env.AppendUnique(LIBPATH = env['CGAL_DIR'])
test_env.AppendUnique(CPPPATH = env['CGAL_INC'])
conf = Configure(test_env)
if conf.CheckLibWithHeader(env['CGAL_LIB'], 'flann/flann.hpp', 'CXX'):
build_lib = True
env.AppendUnique(LIBS = env['CGAL_LIB'])
env.AppendUnique(LIBPATH = env['CGAL_DIR'])
env.AppendUnique(CPPPATH = env['CGAL_INC'])
#end
env = conf.Finish()
if build_lib:
sources = ['data.cpp']
env.SharedLibrary('../../build/data_interpolant', sources)
#end
/* ALTA --- Analysis of Bidirectional Reflectance Distribution Functions
Copyright (C) 2013, 2014 Inria
This file is part of ALTA.
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/. */
#include "data.h"
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cassert>
#include <core/vertical_segment.h>
#ifdef USE_DELAUNAY
#include <CGAL/Cartesian_d.h>
#include <CGAL/Homogeneous_d.h>
#include <CGAL/Delaunay_d.h>
typedef CGAL::Homogeneous_d<double> Kernel;
typedef CGAL::Delaunay_d<Kernel> Delaunay_d;
typedef Delaunay_d::Point_d Point;
typedef Delaunay_d::Simplex_handle S_handle;
typedef Delaunay_d::Point_const_iterator P_iter;
Delaunay_d* D;
int dD;
#else
#include <flann/flann.hpp>
#endif
data_interpolant::data_interpolant() : _data(new vertical_segment())
{
_knn = 3;
}
data_interpolant::~data_interpolant()
{
#ifndef USE_DELAUNAY
if(_kdtree != NULL)
delete _kdtree;
#else
if(D != NULL)
delete D;
#endif
}
// Load data from a file
void data_interpolant::load(const std::string& filename)
{
// Load the data
_data->load(filename);
// Copy the informations
setDimX(_data->dimX());
setDimY(_data->dimY());
setMin(_data->min());
setMax(_data->max());
setParametrization(_data->input_parametrization());
setParametrization(_data->output_parametrization());
#ifdef USE_DELAUNAY
dD = dimX()+dimY();
D = new Delaunay_d(dD);
for(int i=0; i<_data->size(); ++i)
{
vec x = _data->get(i);
Point pt(dD, &x[0], &x[dD]);
D->insert(pt);
}
std::cout << "<<DEBUG>> number of points in the Delaunay triangulation: " << D->all_points().size() << std::endl;
std::cout << "<<DEBUG>> number of points in input: " << _data->size() << std::endl;
#else
// Update the KDtreee by inserting all points
double* _d = new double[dimX()*_data->size()];
flann::Matrix<double> pts(_d, _data->size(), dimX());
for(int i=0; i<_data->size(); ++i)
{
vec x = _data->get(i);
memcpy(pts[i], &x[0], dimX()*sizeof(double));
}
_kdtree = new flann::Index< flann::L2<double> >(pts, flann::KDTreeIndexParams(4));
_kdtree->buildIndex();
#endif
}
void data_interpolant::load(const std::string& filename, const arguments&)
{
load(filename);
}
void data_interpolant::save(const std::string& filename) const
{
}
// Acces to data
vec data_interpolant::get(int id) const
{
vec res(dimX() + dimY()) ;
return res ;
}
vec data_interpolant::operator[](int i) const
{
return get(i) ;
}
//! \todo Test this function
void data_interpolant::set(const vec& x)
{
NOT_IMPLEMENTED();
}
void data_interpolant::set(int i, const vec& x)
{
NOT_IMPLEMENTED();
}
vec data_interpolant::value(const vec& x) const
{
vec res = vec::Zero(dimY());
#ifndef USE_DELAUNAY
// Query point
vec xc(x);
flann::Matrix<double> pts(&xc[0], 1, dimX());
std::vector< std::vector<int> > indices;
std::vector< std::vector<double> > dists;
_kdtree->knnSearch(pts, indices, dists, _knn, flann::SearchParams());
// Interpolate the value using the indices
double cum_dist = 0.0;
for(int i=0; i<indices[0].size(); ++i)
{
const int indice = indices[0][i];
vec y = _data->get(indice);
const double kernel = 1.0/(1.0E-10 + dists[0][i]);
res += kernel * y.tail(_nY);
cum_dist += kernel;
}
if(cum_dist > 0.0)
{
res /= cum_dist;
}
#else
Point pt_x(dD);
for(int j=0; j<dimX(); ++j) { pt_x[j] = x[j]; }
S_handle simplex = D->locate(pt_x);
if(simplex == Delaunay_d::Simplex_handle())
{
return res;
}
// Interpolate the value of the vertices of the Simplex to estimate
// the value of x
double cum_dist = 0.0;
for(int j=0; j<=dD; ++j)
{
Point y = D->point_of_simplex(simplex, j);
// Compute the distance between y and x
double dist = 0.0;
for(int i=0; i<dimX(); ++i) { dist += pow(y.homogeneous(i)-x[i], 2); }
dist = sqrt(dist);
// Interpolate the data
cum_dist += dist;
for(int i=0; i<dimY(); ++i)
{
res[i] += dist * y.homogeneous(dimX() + i);
}
}
if(cum_dist > 0.0)
{
for(int j=0; j<dimY(); ++j)
{
res[j] /= cum_dist;
}
}
#endif
return res;
}
// Get data size, e.g. the number of samples to fit
int data_interpolant::size() const
{
assert(_data);
return _data->size();
}
ALTA_DLL_EXPORT data* provide_data(const arguments&)
{
return new data_interpolant();
}
/* ALTA --- Analysis of Bidirectional Reflectance Distribution Functions
Copyright (C) 2013, 2014 Inria
This file is part of ALTA.
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/. */
#pragma once
#include <core/data.h>
#include <core/common.h>
#include <core/args.h>
//#define USE_DELAUNAY
#ifndef USE_DELAUNAY
#include <flann/flann.hpp>
#endif
/*! \brief Load a data file, but provide access to an interpolated version of
* the data points.
*
* <h2>Requirements:</h2>
* The FLANN library.
* On linux plateforms it can be obtained using the package manager:
* \verbatim
* sudo apt-get install libflann-dev
* \endverbatim
*/
class data_interpolant : public data
{
public: // methods
data_interpolant();
~data_interpolant();
// Load data from a file
virtual void load(const std::string& filename) ;
virtual void load(const std::string& filename, const arguments& args) ;
virtual void save(const std::string& filename) const ;
// Acces to data
virtual vec get(int i) const ;
virtual vec operator[](int i) const ;
virtual vec value(const vec& x) const;
// Set data
virtual void set(const vec& x);
virtual void set(int i, const vec& x);
// Get data size, e.g. the number of samples to fit
virtual int size() const ;
private: // data
// The data object used to load sparse points sets
ptr<data> _data;
// Interpolation
#ifndef USE_DELAUNAY
flann::Index< flann::L2<double> >* _kdtree;
#endif
int _knn;
} ;
/* ALTA --- Analysis of Bidirectional Reflectance Distribution Functions
Copyright (C) 2013 Inria
This file is part of ALTA.
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/. */
#include "data.h"
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cassert>
#include <core/vertical_segment.h>
mxArray *X, *Y, *x, *y;
#define BUFFER_SIZE 10000
char* output = new char[BUFFER_SIZE+1];
data_interpolant::data_interpolant()
{
_data = new vertical_segment();
// Create matlab engine
#ifdef WIN32
if (!(ep = engOpen(NULL)))
#else
if (!(ep = engOpen("matlab -nosplash")))
#endif
{
std::cerr << "<ERROR>> can't start MATLAB engine" << std::endl ;
}
output[BUFFER_SIZE] = '\0';
engOutputBuffer(ep, output, BUFFER_SIZE) ;
}
data_interpolant::~data_interpolant()
{
delete _data;
delete[] output;
mxDestroyArray(x);
mxDestroyArray(y);
mxDestroyArray(X);
mxDestroyArray(Y);
engClose(ep);
}
// Load data from a file
void data_interpolant::load(const std::string& filename)
{
// Load the data
_data->load(filename);
// Copy the informations
setDimX(_data->dimX());
setDimY(_data->dimY());
setMin(_data->min());
setMax(_data->max());
Eigen::MatrixXd eX(_data->size(), dimX()*dimY());
Eigen::MatrixXd eY(_data->size(), dimY());
for(int i=0; i<_data->size(); ++i)
{
vec x = _data->get(i);
// For a point extract all its coordinates
for(int j=0; j<dimX(); ++j)
{
eX(i, j) = x[j];
}
// Extract all its values
for(int k=0; k<dimY(); ++k)
{
eY(i, k) = x[dimX() + k];
}
}
// Create matrices
X = mxCreateDoubleMatrix(_data->size(), dimX(), mxREAL);
memcpy((void *)mxGetPr(X), (void *) eX.data(), _data->size()*dimY()*dimX()*sizeof(double));
engPutVariable(ep, "X", X);
Y = mxCreateDoubleMatrix(_data->size(), dimY(), mxREAL);
memcpy((void *)mxGetPr(Y), (void *) eY.data(), _data->size()*dimY()*sizeof(double));
engPutVariable(ep, "Y", Y);
x = mxCreateDoubleMatrix(1, dimX(), mxREAL);
}
void data_interpolant::load(const std::string& filename, const arguments&)
{
load(filename);
}
void data_interpolant::save(const std::string& filename) const
{
}
// Acces to data
vec data_interpolant::get(int id) const
{
vec res(dimX() + dimY()) ;
return res ;
}
vec data_interpolant::operator[](int i) const
{
return get(i) ;
}
void data_interpolant::set(const vec& x)
{
NOT_IMPLEMENTED();
}
void data_interpolant::set(int i, const vec& x)
{
NOT_IMPLEMENTED();
}
vec data_interpolant::value(const vec& ax) const
{
vec res = vec::Zero(dimY());
// Copy the input vector to matlab code
memcpy((void *)mxGetPr(x), (void *)&ax[0], dimX()*sizeof(double));
engPutVariable(ep, "x", x);
std::stringstream cmd;
// Iterate over the output dimension
for(int i=0; i<dimY(); ++i)
{
// Evaluate the matlab routine
if(dimX() == 2)
{
cmd << "y(" << i+1 << ") = griddata(X(:,1), X(:,2), Y(:," << i+1 << "), x(1), x(2), 'cubic');";
}
else if(dimX() == 3)
{
cmd << "y(" << i+1 << ") = griddata(X(:,1), X(:,2), X(:,3), Y(:," << i+1 << "), x(1), x(2), x(3), 'cubic');";
}
else
{
cmd << "y(" << i+1 << ") = griddatan(X(:, 1:" << dimX() <<"), Y(:, " << i+1 << "), x, 'linear');";
}
engEvalString(ep, cmd.str().c_str());
#ifdef DEBUG
std::cout << output;
#endif
}
// Get results and copy it
y = engGetVariable(ep, "y") ;
double* y_val = (double*)mxGetData(y);
memcpy(&res[0], y_val, dimY()*sizeof(double));
// Fail safe: if the query point is outside of the convex hull of the
// data, rerun the algorithm using a nearest option.
if(isnan(res[0]))
{
for(int i=0; i<dimY(); ++i)
{
cmd.str(std::string());
cmd << "y("<< i+1 <<") = griddatan(X(:, 1:" << dimX() <<"), Y(:, " << i+1 << "), x, 'nearest');";
engEvalString(ep, cmd.str().c_str());
}
// Get results and copy it
y = engGetVariable(ep, "y") ;
double* y_val = (double*)mxGetData(y);
memcpy(&res[0], y_val, dimY()*sizeof(double));
}
return res;
}
// Get data size, e.g. the number of samples to fit
int data_interpolant::size() const
{
assert(_data != NULL);
return _data->size();
}
ALTA_DLL_EXPORT data* provide_data(const arguments&)
{
return new data_interpolant();
}
/* ALTA --- Analysis of Bidirectional Reflectance Distribution Functions
Copyright (C) 2013, 2014 Inria
This file is part of ALTA.
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/. */
#pragma once
#include <core/data.h>
#include <core/common.h>
#include <core/args.h>
#include <engine.h>
/*! \brief Load a data file, but provide access to an interpolated version of
* the data points.
*
* <h2>Requirements:</h2>
* Matlab engine
*/
class data_interpolant : public data
{
public: // methods
data_interpolant();
~data_interpolant();
// Load data from a file
virtual void load(const std::string& filename) ;
virtual void load(const std::string& filename, const arguments& args) ;
virtual void save(const std::string& filename) const ;
// Acces to data
virtual vec get(int i) const ;
virtual vec operator[](int i) const ;
virtual vec value(const vec& x) const;
// Set data
virtual void set(const vec& x);
virtual void set(int i, const vec& x);
// Get data size, e.g. the number of samples to fit
virtual int size() const ;
private: // data
// The data object used to load sparse points sets
ptr<data> _data;
Engine *ep;
} ;
......@@ -4,4 +4,47 @@ env = env.Clone()
# Special linking flags, defined in the OS dependant configuration file
env.AppendUnique(LIBS = env['PLUGIN_LIB'])
build_rbf_lib = False
build_mat_lib = False
test_env = env.Clone()
test_env.AppendUnique(LIBS = env['FLANN_LIB'])
test_env.AppendUnique(LIBPATH = env['FLANN_DIR'])
test_env.AppendUnique(CPPPATH = env['FLANN_INC'])
test_env.AppendUnique(LIBS = env['CGAL_LIB'])
test_env.AppendUnique(LIBPATH = env['CGAL_DIR'])
test_env.AppendUnique(CPPPATH = env['CGAL_INC'])
test_env.AppendUnique(LIBS = env['MATLAB_LIB'])
test_env.AppendUnique(LIBPATH = env['MATLAB_DIR'])
test_env.AppendUnique(CPPPATH = env['MATLAB_INC'])
conf = Configure(test_env)
# Check for the presence of FLANN lib
if conf.CheckLibWithHeader(env['FLANN_LIB'], 'flann/flann.hpp', 'CXX'):
build_rbf_lib = True
env.AppendUnique(LIBS = env['FLANN_LIB'])
env.AppendUnique(LIBPATH = env['FLANN_DIR'])
env.AppendUnique(CPPPATH = env['FLANN_INC'])
# Check for the presence of CGAL lib
if conf.CheckLibWithHeader(env['CGAL_LIB'], 'flann/flann.hpp', 'CXX'):
build_rbf_lib = True
env.AppendUnique(LIBS = env['CGAL_LIB'])
env.AppendUnique(LIBPATH = env['CGAL_DIR'])
env.AppendUnique(CPPPATH = env['CGAL_INC'])
# Check for the presence of Matlab lib
if conf.CheckLibWithHeader(env['MATLAB_LIB'], 'engine.h', 'CXX'):
build_mat_lib = True
env.AppendUnique(LIBS = env['MATLAB_LIB'])
env.AppendUnique(LIBPATH = env['MATLAB_DIR'])
env.AppendUnique(CPPPATH = env['MATLAB_INC'])
test_env = conf.Finish()
env.SharedLibrary('../../build/data_grid', ['grid.cpp'])
if(build_rbf_lib):
env.SharedLibrary('../../build/data_interpolant_rbf', ['rbf.cpp'])
if(build_mat_lib):
env.SharedLibrary('../../build/data_interpolant_matlab', ['matlab.cpp'])
\ No newline at end of file
/* ALTA --- Analysis of Bidirectional Reflectance Distribution Functions
Copyright (C) 2013 Inria
This file is part of ALTA.
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/. */
// STL includes
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cassert>
// ALTA includes
#include <core/data.h>
#include <core/vertical_segment.h>
#include <core/common.h>
#include <core/args.h>
// Matlab includes
#include <engine.h>
mxArray *X, *Y, *x, *y;
#define BUFFER_SIZE 10000
char* output = new char[BUFFER_SIZE+1];
/*! \ingroup datas
* \ingroup plugins
* \class data_interpolant_matlab
* \brief This plugin provide an interpolation of \ref vertical_segment
* data using Matlab `griddata` functions.
*
* \details
* This plugin provide an interpolation of \ref vertical_segment
* data using Matlab `griddata` functions.
*
* ### Requirements
* This plugin requires the Matlab engine library to compile.
*
* \author Laurent Belcour <laurent.belcour@umontreal.ca>
*/
class MatlabInterpolant : public data
{
private: // data
// The data object used to load sparse points sets
ptr<data> _data;
Engine *ep;
public: // methods
MatlabInterpolant()
{
_data = new vertical_segment();
// Create matlab engine
#ifdef WIN32
if (!(ep = engOpen(NULL)))
#else
if (!(ep = engOpen("matlab -nosplash")))
#endif
{
std::cerr << "<ERROR>> can't start MATLAB engine" << std::endl ;
}
output[BUFFER_SIZE] = '\0';
engOutputBuffer(ep, output, BUFFER_SIZE) ;
}
virtual MatlabInterpolant()
{
delete _data;
delete[] output;
mxDestroyArray(x);
mxDestroyArray(y);
mxDestroyArray(X);
mxDestroyArray(Y);
engClose(ep);
}
// Load data from a file
virtual void load(const std::string& filename)
{