rbf.cpp 5.81 KB
Newer Older
1 2
/* ALTA --- Analysis of Bidirectional Reflectance Distribution Functions

3
   Copyright (C) 2013, 2014, 2016 Inria
4 5 6 7 8 9 10 11 12 13 14

   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 <core/data.h>
#include <core/vertical_segment.h>
#include <core/common.h>
#include <core/args.h>
15
#include <core/plugins_manager.h>
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

//#define USE_DELAUNAY
#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

36 37
using namespace alta;

38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
/*! \ingroup datas
 *  \ingroup plugins
 *  \class data_interpolant_rbf
 *  \brief This plugin provide an interpolation of \ref vertical_segment
 *  data using *Radial Basis Function* interpolation.
 *
 *  \details
 *  This interpolation plugin smooth sparse data measurments using radial
 *  kernels. The interpolation is performed as:
 *  <center>
 *     \f$ \tilde{y}(\mathbf{x}) = {\sum_{\mathbf{x}_i \in \mathcal{B}(
 *     \mathbf{x})} k(\mathbf{x} - \mathbf{x}_i) y_i \over \sum_{\mathbf{x}_i
 *     \in \mathcal{B}(\mathbf{x})} k(\mathbf{x} - \mathbf{x}_i)}\f$
 *  </center>
 *
 *  The kernel is set as \f$ k(\mathbf{d}) = {1 \over \epsilon + ||d||^2} \f$.
 *
 *  ### Requirements
 *  This plugin requires the FLANN library to compile. On linux plateforms
 *  it can be obtained by package `libflann-dev` and on OSX using the port
 *  `flann`.
 *
 *  \author Laurent Belcour <laurent.belcour@umontreal.ca>
 */
class rbf_interpolant : public data
{
	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;

	public:

76 77 78 79
		rbf_interpolant(ptr<data> proxied_data)
        : data(proxied_data->parametrization()),
          _data(proxied_data),
          _knn(3)
80
		{
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
        _min = _data->min();
        _max = _data->max();

#ifdef USE_DELAUNAY
        dD = parametrization().dimX()+parametrization().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[parametrization().dimX()*_data->size()];
        flann::Matrix<double> pts(_d, _data->size(), parametrization().dimX());
        for(int i=0; i<_data->size(); ++i)
        {
            vec x = _data->get(i);
            memcpy(pts[i], &x[0], parametrization().dimX()*sizeof(double));
        }
        _kdtree = new flann::Index< flann::L2<double> >(pts, flann::KDTreeIndexParams(4));
        _kdtree->buildIndex();
#endif
    }
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128

		virtual ~rbf_interpolant()
		{
		#ifndef USE_DELAUNAY
			if(_kdtree != NULL)
				delete _kdtree;
		#else
			if(D != NULL)
				delete D;
		#endif
		}

		virtual void save(const std::string& filename) const
		{
		}

		// Acces to data
		virtual vec get(int id) const
		{
129
			vec res(parametrization().dimX() + parametrization().dimY()) ;
130 131 132 133 134 135 136 137 138 139 140 141 142 143
			return res ;
		}
		virtual vec operator[](int i) const
		{
			return get(i) ;
		}

		virtual void set(int i, const vec& x)
		{
			NOT_IMPLEMENTED();
		}

		virtual vec value(const vec& x) const
		{
144
			vec res = vec::Zero(parametrization().dimY());
145 146 147 148

		#ifndef USE_DELAUNAY
			// Query point
			vec xc(x);
149
			flann::Matrix<double> pts(&xc[0], 1, parametrization().dimX());
150 151 152 153 154 155 156 157 158 159 160 161 162 163
			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]);

164
        res      += kernel * y.tail(parametrization().dimY());
165 166 167 168 169 170 171 172 173
				cum_dist += kernel;
			}
			if(cum_dist > 0.0)
			{
				res /= cum_dist;
			}
		#else

			Point pt_x(dD);
174
			for(int j=0; j<parametrization().dimX(); ++j) { pt_x[j] = x[j]; }
175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
			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;
191
				for(int i=0; i<parametrization().dimX(); ++i) { dist += pow(y.homogeneous(i)-x[i], 2); }
192 193 194 195
				dist = sqrt(dist);

				// Interpolate the data
				cum_dist += dist;
196
				for(int i=0; i<parametrization().dimY(); ++i)
197
				{
198
					res[i] += dist * y.homogeneous(parametrization().dimX() + i);
199 200 201 202 203 204
				}

			}

			if(cum_dist > 0.0)
			{
205
				for(int j=0; j<parametrization().dimY(); ++j)
206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222
				{
					res[j] /= cum_dist;
				}
			}
		#endif

		   return res;
		}

		// Get data size, e.g. the number of samples to fit
		virtual int size() const
		{
			assert(_data);
			return _data->size();
		}
};

223
ALTA_DLL_EXPORT data* load_data(std::istream& input, const arguments& args)
224
{
225 226 227 228 229
    // Load the data
    ptr<data> proxied = plugins_manager::load_data("vertical_segment",
                                                   input, args);

    return new rbf_interpolant(proxied);
230
}
231 232