有没有什么好的库可以在 C++ 中解决三次样条？

Are there any good libraries for solving cubic splines in C++?

提问by Faken

I'm looking for a good C++ library to give me functions to solve for large cubic splines (on the order of 1000 points) anyone know one?

我正在寻找一个好的 C++ 库来为我提供解决大型三次样条（大约 1000 点）的函数，有人知道吗？

采纳答案by ars

Try the Cubic B-Spline library:

试试 Cubic B-Spline 库：

• https://github.com/NCAR/bspline

• https://github.com/NCAR/bspline

and ALGLIB:

和 ALGLIB：

• http://www.alglib.net/interpolation/spline3.php

• http://www.alglib.net/interpolation/spline3.php

回答by cpp

Write your own. Here is spline()function I wrote based on excellent wiki algorithm:

自己写。这是spline()我根据优秀的维基算法编写的函数：

#include<iostream>
#include<vector>
#include<algorithm>
#include<cmath>
using namespace std;

using vec = vector<double>;

struct SplineSet{
    double a;
    double b;
    double c;
    double d;
    double x;
};

vector<SplineSet> spline(vec &x, vec &y)
{
    int n = x.size()-1;
    vec a;
    a.insert(a.begin(), y.begin(), y.end());
    vec b(n);
    vec d(n);
    vec h;

    for(int i = 0; i < n; ++i)
        h.push_back(x[i+1]-x[i]);

    vec alpha;
    alpha.push_back(0);
    for(int i = 1; i < n; ++i)
        alpha.push_back( 3*(a[i+1]-a[i])/h[i] - 3*(a[i]-a[i-1])/h[i-1]  );

    vec c(n+1);
    vec l(n+1);
    vec mu(n+1);
    vec z(n+1);
    l[0] = 1;
    mu[0] = 0;
    z[0] = 0;

    for(int i = 1; i < n; ++i)
    {
        l[i] = 2 *(x[i+1]-x[i-1])-h[i-1]*mu[i-1];
        mu[i] = h[i]/l[i];
        z[i] = (alpha[i]-h[i-1]*z[i-1])/l[i];
    }

    l[n] = 1;
    z[n] = 0;
    c[n] = 0;

    for(int j = n-1; j >= 0; --j)
    {
        c[j] = z [j] - mu[j] * c[j+1];
        b[j] = (a[j+1]-a[j])/h[j]-h[j]*(c[j+1]+2*c[j])/3;
        d[j] = (c[j+1]-c[j])/3/h[j];
    }

    vector<SplineSet> output_set(n);
    for(int i = 0; i < n; ++i)
    {
        output_set[i].a = a[i];
        output_set[i].b = b[i];
        output_set[i].c = c[i];
        output_set[i].d = d[i];
        output_set[i].x = x[i];
    }
    return output_set;
}

int main()
{
    vec x(11);
    vec y(11);
    for(int i = 0; i < x.size(); ++i)
    {
        x[i] = i;
        y[i] = sin(i);
    }

    vector<SplineSet> cs = spline(x, y);
    for(int i = 0; i < cs.size(); ++i)
        cout << cs[i].d << "\t" << cs[i].c << "\t" << cs[i].b << "\t" << cs[i].a << endl;
}

回答by FuzzyBunnySlippers

I had to write spline routine for an "entity" that was following a path (series of connected waypoints) in a game I am working on.

我必须为一个“实体”编写样条例程，该实体在我正在开发的游戏中遵循路径（一系列连接的航点）。

I created a base class to handle a "SplineInterface" and the created two derived classes, one based on the classic spline technique (e.g. Sedgewick/Algorithms) an a second one based on Bezier Splines.

我创建了一个基类来处理“SplineInterface”，并创建了两个派生类，一个基于经典样条技术（例如 Sedgewick/Algorithms），另一个基于 Bezier Splines。

Here is the code. It is a single header file, which contains all the splining classes:

这是代码。它是一个单独的头文件，其中包含所有样条类：

#ifndef __SplineCommon__
#define __SplineCommon__

#include "CommonSTL.h"
#include "CommonProject.h"
#include "MathUtilities.h"

/* A Spline base class. */
class SplineBase
{
private:
   vector<Vec2> _points;
   bool _elimColinearPoints;

protected:


protected:
   /* OVERRIDE THESE FUNCTIONS */
   virtual void ResetDerived() = 0;

   enum
   {
      NOM_SIZE = 32,
   };

public:

   SplineBase()
   {
      _points.reserve(NOM_SIZE);
      _elimColinearPoints = true;
   }

   const vector<Vec2>& GetPoints() { return _points; }
   bool GetElimColinearPoints() { return _elimColinearPoints; }
   void SetElimColinearPoints(bool elim) { _elimColinearPoints = elim; }


   /* OVERRIDE THESE FUNCTIONS */
   virtual Vec2 Eval(int seg, double t) = 0;
   virtual bool ComputeSpline() = 0;
   virtual void DumpDerived() {}

   /* Clear out all the data.
    */
   void Reset()
   {
      _points.clear();
      ResetDerived();
   }

   void AddPoint(const Vec2& pt)
   {
      // If this new point is colinear with the two previous points,
      // pop off the last point and add this one instead.
      if(_elimColinearPoints && _points.size() > 2)
      {
         int N = _points.size()-1;
         Vec2 p0 = _points[N-1] - _points[N-2];
         Vec2 p1 = _points[N] - _points[N-1];
         Vec2 p2 = pt - _points[N];
         // We test for colinearity by comparing the slopes
         // of the two lines.  If the slopes are the same,
         // we assume colinearity.
         float32 delta = (p2.y-p1.y)*(p1.x-p0.x)-(p1.y-p0.y)*(p2.x-p1.x);
         if(MathUtilities::IsNearZero(delta))
         {
            _points.pop_back();
         }
      }
      _points.push_back(pt);
   }

   void Dump(int segments = 5)
   {
      assert(segments > 1);

      cout << "Original Points (" << _points.size() << ")" << endl;
      cout << "-----------------------------" << endl;
      for(int idx = 0; idx < _points.size(); ++idx)
      {
         cout << "[" << idx << "]" << "  " << _points[idx] << endl;
      }

      cout << "-----------------------------" << endl;
      DumpDerived();

      cout << "-----------------------------" << endl;
      cout << "Evaluating Spline at " << segments << " points." << endl;
      for(int idx = 0; idx < _points.size()-1; idx++)
      {
         cout << "---------- " << "From " <<  _points[idx] << " to " << _points[idx+1] << "." << endl;
         for(int tIdx = 0; tIdx < segments+1; ++tIdx)
         {
            double t = tIdx*1.0/segments;
            cout << "[" << tIdx << "]" << "   ";
            cout << "[" << t*100 << "%]" << "   ";
            cout << " --> " << Eval(idx,t);
            cout << endl;
         }
      }
   }
};

class ClassicSpline : public SplineBase
{
private:
   /* The system of linear equations found by solving
    * for the 3 order spline polynomial is given by:
    * A*x = b.  The "x" is represented by _xCol and the
    * "b" is represented by _bCol in the code.
    *
    * The "A" is formulated with diagonal elements (_diagElems) and
    * symmetric off-diagonal elements (_offDiagElemns).  The
    * general structure (for six points) looks like:
    *
    *
    *  |  d1  u1   0   0   0  |      | p1 |    | w1 |
    *  |  u1  d2   u2  0   0  |      | p2 |    | w2 |
    *  |  0   u2   d3  u3  0  |   *  | p3 |  = | w3 |
    *  |  0   0    u3  d4  u4 |      | p4 |    | w4 |
    *  |  0   0    0   u4  d5 |      | p5 |    | w5 |
    *
    *
    *  The general derivation for this can be found
    *  in Robert Sedgewick's "Algorithms in C++".
    *
    */
   vector<double> _xCol;
   vector<double> _bCol;
   vector<double> _diagElems;
   vector<double> _offDiagElems;
public:
   ClassicSpline()
   {
      _xCol.reserve(NOM_SIZE);
      _bCol.reserve(NOM_SIZE);
      _diagElems.reserve(NOM_SIZE);
      _offDiagElems.reserve(NOM_SIZE);
   }

   /* Evaluate the spline for the ith segment
    * for parameter.  The value of parameter t must
    * be between 0 and 1.
    */
   inline virtual Vec2 Eval(int seg, double t)
   {
      const vector<Vec2>& points = GetPoints();

      assert(t >= 0);
      assert(t <= 1.0);
      assert(seg >= 0);
      assert(seg < (points.size()-1));

      const double ONE_OVER_SIX = 1.0/6.0;
      double oneMinust = 1.0 - t;
      double t3Minust = t*t*t-t;
      double oneMinust3minust = oneMinust*oneMinust*oneMinust-oneMinust;
      double deltaX = points[seg+1].x - points[seg].x;
      double yValue = t * points[seg + 1].y +
      oneMinust*points[seg].y +
      ONE_OVER_SIX*deltaX*deltaX*(t3Minust*_xCol[seg+1] - oneMinust3minust*_xCol[seg]);
      double xValue = t*(points[seg+1].x-points[seg].x) + points[seg].x;
      return Vec2(xValue,yValue);
   }


   /* Clear out all the data.
    */
   virtual void ResetDerived()
   {
      _diagElems.clear();
      _bCol.clear();
      _xCol.clear();
      _offDiagElems.clear();
   }


   virtual bool ComputeSpline()
   {
      const vector<Vec2>& p = GetPoints();


      _bCol.resize(p.size());
      _xCol.resize(p.size());
      _diagElems.resize(p.size());

      for(int idx = 1; idx < p.size(); ++idx)
      {
         _diagElems[idx] = 2*(p[idx+1].x-p[idx-1].x);
      }
      for(int idx = 0; idx < p.size(); ++idx)
      {
         _offDiagElems[idx] = p[idx+1].x - p[idx].x;
      }
      for(int idx = 1; idx < p.size(); ++idx)
      {
         _bCol[idx] = 6.0*((p[idx+1].y-p[idx].y)/_offDiagElems[idx] -
                           (p[idx].y-p[idx-1].y)/_offDiagElems[idx-1]);
      }
      _xCol[0] = 0.0;
      _xCol[p.size()-1] = 0.0;
      for(int idx = 1; idx < p.size()-1; ++idx)
      {
         _bCol[idx+1] = _bCol[idx+1] - _bCol[idx]*_offDiagElems[idx]/_diagElems[idx];
         _diagElems[idx+1] = _diagElems[idx+1] - _offDiagElems[idx]*_offDiagElems[idx]/_diagElems[idx];
      }
      for(int idx = (int)p.size()-2; idx > 0; --idx)
      {
         _xCol[idx] = (_bCol[idx] - _offDiagElems[idx]*_xCol[idx+1])/_diagElems[idx];
      }
      return true;
   }
};

/* Bezier Spline Implementation
 * Based on this article:
 * http://www.particleincell.com/blog/2012/bezier-splines/
 */
class BezierSpine : public SplineBase
{
private:
   vector<Vec2> _p1Points;
   vector<Vec2> _p2Points;
public:
   BezierSpine()
   {
      _p1Points.reserve(NOM_SIZE);
      _p2Points.reserve(NOM_SIZE);
   }

   /* Evaluate the spline for the ith segment
    * for parameter.  The value of parameter t must
    * be between 0 and 1.
    */
   inline virtual Vec2 Eval(int seg, double t)
   {
      assert(seg < _p1Points.size());
      assert(seg < _p2Points.size());

      double omt = 1.0 - t;

      Vec2 p0 = GetPoints()[seg];
      Vec2 p1 = _p1Points[seg];
      Vec2 p2 = _p2Points[seg];
      Vec2 p3 = GetPoints()[seg+1];

      double xVal = omt*omt*omt*p0.x + 3*omt*omt*t*p1.x +3*omt*t*t*p2.x+t*t*t*p3.x;
      double yVal = omt*omt*omt*p0.y + 3*omt*omt*t*p1.y +3*omt*t*t*p2.y+t*t*t*p3.y;
      return Vec2(xVal,yVal);
   }

   /* Clear out all the data.
    */
   virtual void ResetDerived()
   {
      _p1Points.clear();
      _p2Points.clear();
   }


   virtual bool ComputeSpline()
   {
      const vector<Vec2>& p = GetPoints();

      int N = (int)p.size()-1;
      _p1Points.resize(N);
      _p2Points.resize(N);
      if(N == 0)
         return false;

      if(N == 1)
      {  // Only 2 points...just create a straight line.
         // Constraint:  3*P1 = 2*P0 + P3
         _p1Points[0] = (2.0/3.0*p[0] + 1.0/3.0*p[1]);
         // Constraint:  P2 = 2*P1 - P0
         _p2Points[0] = 2.0*_p1Points[0] - p[0];
         return true;
      }

      /*rhs vector*/
      vector<Vec2> a(N);
      vector<Vec2> b(N);
      vector<Vec2> c(N);
      vector<Vec2> r(N);

      /*left most segment*/
      a[0].x = 0;
      b[0].x = 2;
      c[0].x = 1;
      r[0].x = p[0].x+2*p[1].x;

      a[0].y = 0;
      b[0].y = 2;
      c[0].y = 1;
      r[0].y = p[0].y+2*p[1].y;

      /*internal segments*/
      for (int i = 1; i < N - 1; i++)
      {
         a[i].x=1;
         b[i].x=4;
         c[i].x=1;
         r[i].x = 4 * p[i].x + 2 * p[i+1].x;

         a[i].y=1;
         b[i].y=4;
         c[i].y=1;
         r[i].y = 4 * p[i].y + 2 * p[i+1].y;
      }

      /*right segment*/
      a[N-1].x = 2;
      b[N-1].x = 7;
      c[N-1].x = 0;
      r[N-1].x = 8*p[N-1].x+p[N].x;

      a[N-1].y = 2;
      b[N-1].y = 7;
      c[N-1].y = 0;
      r[N-1].y = 8*p[N-1].y+p[N].y;


      /*solves Ax=b with the Thomas algorithm (from Wikipedia)*/
      for (int i = 1; i < N; i++)
      {
         double m;

         m = a[i].x/b[i-1].x;
         b[i].x = b[i].x - m * c[i - 1].x;
         r[i].x = r[i].x - m * r[i-1].x;

         m = a[i].y/b[i-1].y;
         b[i].y = b[i].y - m * c[i - 1].y;
         r[i].y = r[i].y - m * r[i-1].y;
      }

      _p1Points[N-1].x = r[N-1].x/b[N-1].x;
      _p1Points[N-1].y = r[N-1].y/b[N-1].y;
      for (int i = N - 2; i >= 0; --i)
      {
         _p1Points[i].x = (r[i].x - c[i].x * _p1Points[i+1].x) / b[i].x;
         _p1Points[i].y = (r[i].y - c[i].y * _p1Points[i+1].y) / b[i].y;
      }

      /*we have p1, now compute p2*/
      for (int i=0;i<N-1;i++)
      {
         _p2Points[i].x=2*p[i+1].x-_p1Points[i+1].x;
         _p2Points[i].y=2*p[i+1].y-_p1Points[i+1].y;
      }

      _p2Points[N-1].x = 0.5 * (p[N].x+_p1Points[N-1].x);
      _p2Points[N-1].y = 0.5 * (p[N].y+_p1Points[N-1].y);

      return true;
   }

   virtual void DumpDerived()
   {
      cout << " Control Points " << endl;
      for(int idx = 0; idx < _p1Points.size(); idx++)
      {
         cout << "[" << idx << "]  ";
         cout << "P1: " << _p1Points[idx];
         cout << "   ";
         cout << "P2: " << _p2Points[idx];
         cout << endl;
      }
   }
};


#endif /* defined(__SplineCommon__) */

Some Notes

一些注意事项

• The classic spline will crash if you give it a vertical set of points. That is why I created the Bezier...I have lots of vertical lines/paths to follow. It could be modified to just give a straight line.
• The base class has an option to remove collinear points as you add them. This uses a simple slope comparison of two lines to figure out if they are on the same line. You don't have to do this, but for long paths that are straight lines, it cuts down on cycles. When you do a lot of pathfinding on a regular-spaced graph, you tend to get a lot of continuous segments.

• 如果给它一组垂直的点，经典样条线会崩溃。这就是我创建贝塞尔曲线的原因……我有很多垂直线/路径要遵循。它可以修改为只给出一条直线。
• 基类可以选择在您添加共线点时删除它们。这使用两条线的简单斜率比较来确定它们是否在同一条线上。您不必这样做，但对于直线的长路径，它会减少周期。当你在一个规则间隔的图上进行大量的寻路时，你往往会得到很多连续的段。

Here is an example of using the Bezier Spline:

以下是使用 Bezier 样条曲线的示例：

/* Smooth the points on the path so that turns look
 * more natural.  We'll only smooth the first few 
 * points.  Most of the time, the full path will not
 * be executed anyway...why waste cycles.
 */
void SmoothPath(vector<Vec2>& path, int32 divisions)
{
   const int SMOOTH_POINTS = 6;

   BezierSpine spline;

   if(path.size() < 2)
      return;

   // Cache off the first point.  If the first point is removed,
   // the we occasionally run into problems if the collision detection
   // says the first node is occupied but the splined point is too
   // close, so the FSM "spins" trying to find a sensor cell that is
   // not occupied.
   //   Vec2 firstPoint = path.back();
   //   path.pop_back();
   // Grab the points.
   for(int idx = 0; idx < SMOOTH_POINTS && path.size() > 0; idx++)
   {
      spline.AddPoint(path.back());
      path.pop_back();
   }
   // Smooth them.
   spline.ComputeSpline();
   // Push them back in.
   for(int idx = spline.GetPoints().size()-2; idx >= 0; --idx)
   {
      for(int division = divisions-1; division >= 0; --division)
      {
         double t = division*1.0/divisions;
         path.push_back(spline.Eval(idx, t));
      }
   }
   // Push back in the original first point.
   //   path.push_back(firstPoint);
}

Notes

笔记

• While the whole path could be smoothed, in this application, since the path was changing every so often, it was better to just smooth the first points and then connect it up.
• The points are loaded in "reverse" order into the path vector. This may or may not save cycles (I've slept since then).

• 虽然可以平滑整个路径，但在此应用程序中，由于路径经常变化，最好只平滑第一个点，然后将其连接起来。
• 这些点以“反向”顺序加载到路径向量中。这可能会也可能不会节省周期（从那时起我就睡了）。

This code is part of a much larger code base, but you can download it all on githuband see a blog entry about it here.

此代码是一个更大的代码库的一部分，但您可以在 github 上下载所有代码并在此处查看有关它的博客条目。

You can look at this in action in this video.

您可以在此视频中查看实际操作。

回答by denis

Take a look at David Eberly's GeometricTools.com. I'm just starting, but code and doc are so far of outstanding quality.
(He has books too: Geometric tools for computer graphics, 3D game engine design.)

看看 David Eberly 的GeometricTools.com。我才刚刚开始，但到目前为止，代码和文档的质量非常出色。
（他也有书：用于计算机图形学的几何工具、3D 游戏引擎设计。）