使用导数对特征样条进行插值

Question

我正在尝试使用指定的开始和结束导数来插值三次样条。

为此，我为特征样条创建了一个包装器：

/**
* Cubic Spline.
* @tparam vec_t point type. Example: Eigen::Vector2d or Eigen::Vector3d
*/
template <typename vec_t, int dim, int spline_degree = 3>
class CubicSpline {
    private:
    typedef typename Eigen::Spline<double, dim, spline_degree>::KnotVectorType KnotVectorType;

    Eigen::MatrixXd m_Points;                               // Points.
    Eigen::Spline<double, dim, spline_degree> m_Spline;     // The spline object.
    Eigen::Spline<double, dim, spline_degree> m_SplineDer;  // The spline obtained with derivatives.
    KnotVectorType m_ChordLengths;  // Knot parameters, [0, 1].

    public:
    // Constructor.
    CubicSpline(const std::vector<vec_t>& points) {
        // Map data to matrix.
        m_Points = Eigen::MatrixXd::Map(points[0].data(), dim, points.size());
        // Intepolate spline.
        Eigen::ChordLengths(m_Points, m_ChordLengths);
        m_Spline = Eigen::SplineFitting<decltype(m_Spline)>::Interpolate(m_Points, spline_degree, m_ChordLengths);
    };

    // THIS FUNCTION FAILS!
    void withDerivatives()
            {
                Eigen::MatrixXd derivatives(2, 1);
                derivatives.setZero(); // Derivatives zero for this example.

                Eigen::Vector<int, 2> indices;
                indices << 0, m_Points.cols() - 1;
                
                m_SplineDer =
                    Eigen::SplineFitting<decltype(m_SplineDer)>::InterpolateWithDerivatives(
                        m_Points, derivatives, indices, spline_degree, m_ChordLengths);
            };

    // Value.
    constexpr vec_t valueAt(const double& t) { return m_Spline(t); };

};

// Common instantiations.
template class CubicSpline<Eigen::Vector2d, 2>;
template class CubicSpline<Eigen::Vector3d, 3>;

然后我在主cpp中进行测试：

#include <unsupported/Eigen/Splines>
#include <vector>
int main(){
    // interpolation points.
    std::vector<Eigen::Vector2d> points(5);
    points[0] = Eigen::Vector2d{0.0, 0.0};
    points[1] = Eigen::Vector2d{10.0, 10.0};
    points[2] = Eigen::Vector2d{20.0, 20.0};
    points[3] = Eigen::Vector2d{30.0, 30.0};
    points[4] = Eigen::Vector2d{40.0, 40.0};

    auto spline = CubicSpline<Eigen::Vector2d, 2>(points);
    size_t N = 20;
    for (size_t i = 0; i < N; ++i) {
        auto t = static_cast<double>(i) / N;

        // Evaluate some points.
        Eigen::Vector2d obtained = spline.valueAt(t);
        Eigen::Vector2d expected = t * (points[0] + points[4]);
        for (size_t j = 0; j < 2; ++j) {
            if (abs(obtained[j] - expected[j]) > 1e-3){
                printf("NOT GOOD: %i,%i \n", i, j);
            }
        }
    }

    printf("now fit with derivatives:\n");
    spline.withDerivatives(); // <----THIS FAILS
    printf("DONE!");

    return 0;
}

问题

它可以编译，但我收到一个关于列和行大小的奇怪特征错误。

output.s: /opt/compiler-explorer/libs/eigen/v3.4.0/Eigen/src/Core/Block.h:146: Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>::Block(XprType&, Eigen::Index, Eigen::Index, Eigen::Index, Eigen::Index) [with XprType = Eigen::Block<Eigen::Matrix<double, -1, -1>, 1, -1, false>; int BlockRows = 1; int BlockCols = -1; bool InnerPanel = false; Eigen::Index = long int]: Assertion `startRow >= 0 && blockRows >= 0 && startRow <= xpr.rows() - blockRows && startCol >= 0 && blockCols >= 0 && startCol <= xpr.cols() - blockCols' failed.
Program terminated with signal: SIGSEGV

有人知道如何解决这个问题吗？

一个工作示例

可以在这里找到

Answer 1

断言失败来自于这样一个事实：当 构造样条函数设置方程 对于每个插值/导数条件， 这意味着将行写成矩阵

A

和

b

将使系统

A * x = b

，其中

x

是样条线 参数，它尝试在外部的

A

或

b

中写入 它假定这些矩阵的大小。

在这种情况下，问题的原因很简单 第一行

Eigen::MatrixXd derivatives(2, 1);
derivatives.setZero(); // Derivatives zero for this example.

由于这些是参数样条线，您必须提供

x

和

y

的导数，对于两点中的每一个 您已选择，所以正确的尺寸是：

Eigen::MatrixXd derivatives(2, 2);
derivatives.setZero(); // Derivatives zero for this example.

作为奖励:)，我绘制区间 [0, 1] 中的样条值， 对于简单插值和导数的情况， 一种情况是

x(t)

和

y(t)

的零导数 在上图中，由于点不均匀， 但正在向两端“聚集” 导数被设置为零。