You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
66 lines
1.8 KiB
C++
66 lines
1.8 KiB
C++
#ifndef RBF_INTERPOLATOR_HPP
|
|
#define RBF_INTERPOLATOR_HPP
|
|
|
|
#include <Eigen/Dense>
|
|
#include <vector>
|
|
#include <cmath>
|
|
|
|
class RBFInterpolator {
|
|
public:
|
|
RBFInterpolator(const Eigen::MatrixXd& points, const Eigen::VectorXd& values, double smooth = 0.0)
|
|
: m_points(points), m_values(values), m_smooth(smooth)
|
|
{
|
|
computeWeights();
|
|
}
|
|
|
|
double interpolate(double x, double y) const {
|
|
Eigen::VectorXd rbfVec(m_points.rows());
|
|
for (int i = 0; i < m_points.rows(); ++i) {
|
|
double r = (Eigen::Vector2d(x, y) - m_points.row(i).transpose()).norm();
|
|
rbfVec(i) = thinPlateKernel(r);
|
|
}
|
|
return rbfVec.dot(m_weights.head(m_weights.size() - 3))
|
|
+ m_weights.tail(3).dot(Eigen::Vector3d(1, x, y));
|
|
}
|
|
|
|
private:
|
|
Eigen::MatrixXd m_points;
|
|
Eigen::VectorXd m_values;
|
|
Eigen::VectorXd m_weights;
|
|
double m_smooth;
|
|
|
|
static double thinPlateKernel(double r) {
|
|
if (r == 0) return 0.0;
|
|
return r * r * std::log(r);
|
|
}
|
|
|
|
void computeWeights() {
|
|
int N = m_points.rows();
|
|
Eigen::MatrixXd K = Eigen::MatrixXd::Zero(N + 3, N + 3);
|
|
|
|
for (int i = 0; i < N; ++i) {
|
|
for (int j = 0; j < N; ++j) {
|
|
double r = (m_points.row(i) - m_points.row(j)).norm();
|
|
K(i, j) = thinPlateKernel(r);
|
|
}
|
|
// P¾ØÕó
|
|
K(i, N) = 1.0;
|
|
K(i, N + 1) = m_points(i, 0);
|
|
K(i, N + 2) = m_points(i, 1);
|
|
K(N, i) = 1.0;
|
|
K(N + 1, i) = m_points(i, 0);
|
|
K(N + 2, i) = m_points(i, 1);
|
|
}
|
|
|
|
// ƽ»¬´¦Àí
|
|
K.topLeftCorner(N, N) += m_smooth * Eigen::MatrixXd::Identity(N, N);
|
|
|
|
Eigen::VectorXd rhs = Eigen::VectorXd::Zero(N + 3);
|
|
rhs.head(N) = m_values;
|
|
|
|
m_weights = K.fullPivLu().solve(rhs);
|
|
}
|
|
};
|
|
|
|
#endif // RBF_INTERPOLATOR_HPP
|