#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