最近在写一个程序，其中需要对B样条曲线进行拟合。但是B样条曲线的公式实在复杂，看着就头晕。于是，我将问题进行了简化。一段B样条曲线，可以近似地看成是若干段抛物线构成的，所以，曲线拟合问题就被转换为抛物线拟合问题了。对于抛物线拟合问题，可以使用《计算方法》中的最小二乘法，最后求解线性方程组的地方，用的是高斯消去法。本文用C#实现了这两种算法。
最小二乘法是一种数据优化技术，在已经得到一组数据的情况下，通过最小化误差平方和的办法，找出最接近的函数。关于最小二乘法的详细说明，可以查看维基百科，本文直接把其中的公式拿来使用。
假设有一组实测的坐标点数据，绘制在窗体上，效果如下图所示：
从图上可以看出，这非常接近一条抛物线。我们将抛物线定义为：a + b*x + c*x^2 = f(x)
我们可以很方便地计算出：n、Σx、Σx^2、Σx^3、Σx^4、Σy、Σx*y、Σx^2*y。C#代码如下（变量开头的s表示sigma）：
int n = points.Count; double sx = (double)points.Sum(c => c.X); double sx2 = (double)points.Sum(c => Math.Pow(c.X, 2)); double sx3 = (double)points.Sum(c => Math.Pow(c.X, 3)); double sx4 = (double)points.Sum(c => Math.Pow(c.X, 4)); double sy = (double)points.Sum(c => c.Y); double sxy = (double)points.Sum(c => c.X * c.Y); double sx2y = (double)points.Sum(c => Math.Pow(c.X, 2) * c.Y);
根据上述数据，我们可以构造一个线性方程组，其正规形式如下：
对于线性方程组，可以使用高斯消去法来求解。我实现了一个通用的高斯消去算法（Gauss），将上述参数代入后，得到一组唯一解：-557.884 + 5.15*x -0.008*x^2 = f(x)
拟合效果如下图所示（红色是实测数据，蓝色是拟合的抛物线）：
完整的C#源代码如下（感谢“zhuangzhuang1988”同学的提醒，我对Gauss方法进行了改进，加入了主元选择，以提高计算精度）：
View Code
public partial class Form1 : Form { public Form1() { InitializeComponent(); } private string testData = "151,41;153,46;156,53;158,57;162,66;164,73;168,81;173,94;177,102;181,116;186,125;192,138;196,150;203,162;207,172;217,187;223,198;232,208;237,215;245,224;248,230;254,237;258,240;264,245;273,249;276,250;283,252;287,254;318,254;339,250;347,247;359,242;366,237;378,232;385,227;406,212;420,201;429,190;437,180;446,168;458,146;467,127;473,107;474,100;476,85;478,75;480,57;487,34;487,33"; /// <summary> /// 解析实测数据。 /// </summary> /// <param name="data"></param> /// <returns></returns> private List<Point> ParseTestData(string data) { List<Point> points = new List<Point>(); foreach (var pointString in data.Split(';')) { string[] items = pointString.Split(','); points.Add(new Point(int.Parse(items[0]), int.Parse(items[1]))); } return points; } protected override void OnPaint(PaintEventArgs e) { //绘制实测数据 var points = ParseTestData(testData); foreach (var pt in points) { e.Graphics.DrawEllipse(Pens.Red, pt.X - 2, pt.Y - 2, 5, 5); } //抛物线拟合 var fittingPoints = ParaCurveFitting(points); //绘制拟合后的抛物线 e.Graphics.DrawLines(Pens.Blue, fittingPoints.ToArray()); } /// <summary> /// 抛物线拟合。 /// </summary> /// <param name="points">实测数据。</param> /// <returns>拟合数据。</returns> private List<Point> ParaCurveFitting(List<Point> points) { if (points.Count < 2) return points; //构造线性方程组 //a + bx + cx^2 = f(x) int n = points.Count; double sx = (double)points.Sum(c => c.X); double sx2 = (double)points.Sum(c => Math.Pow(c.X, 2)); double sx3 = (double)points.Sum(c => Math.Pow(c.X, 3)); double sx4 = (double)points.Sum(c => Math.Pow(c.X, 4)); double sy = (double)points.Sum(c => c.Y); double sxy = (double)points.Sum(c => c.X * c.Y); double sx2y = (double)points.Sum(c => Math.Pow(c.X, 2) * c.Y); //高斯消去法求解方程组 var result = Gauss(new double[,] { { n, sx, sx2, sy }, { sx, sx2, sx3, sxy }, { sx2, sx3, sx4, sx2y } }); //拟合后的抛物线数据 List<Point> fittingPoints = new List<Point>(); if (result != null && result.Length >= 2) { double a = result[0]; double b = result[1]; double c = result[2]; foreach (var pt in points) { fittingPoints.Add(new Point(pt.X, (int)(a + b * pt.X + c * pt.X * pt.X))); } } return fittingPoints; } /// <summary> /// 高斯消去法求线性方程组的解。 /// </summary> /// <param name="matrix">增广系数矩阵。</param> /// <returns>方程组的解（从低次到高次）。如果无解，则返回空引用；如果有无穷解，则返回空数组。</returns> /// <exception cref="ArgumentNullException">如果matrix为空引用，则抛出该异常。</exception> static double[] Gauss(double[,] matrix) { if (matrix == null) throw new ArgumentNullException("matrix"); //无解 int cols = matrix.GetLength(1); if (cols <= 1) return null; //有无穷解 int rows = matrix.GetLength(0); if (rows < cols - 1) return new double[] { }; //转换为行阶梯 for (int i = 0; i < rows - 1; i++) { //选取主元（提高计算精度） GaussPivoting(matrix, i, rows, cols); //消去一列 for (int k = i + 1; k < rows; k++) { if (matrix[k, i] != 0) { double t = matrix[i, i] / matrix[k, i]; for (int j = i; j < cols - 1; j++) { matrix[k, j] = matrix[i, j] - matrix[k, j] * t; } matrix[k, cols - 1] = matrix[i, cols - 1] - matrix[k, cols - 1] * t; } } } //检查秩，判断是否有唯一解 int rank1 = 0, rank2 = 0; for (int i = 0; i < rows; i++) { bool isZeroRow = true; for (int j = i; j < cols - 1; j++) { if (matrix[i, j] != 0) { isZeroRow = false; break; } } if (!isZeroRow) rank1++; if (!isZeroRow || matrix[i, cols - 1] != 0) rank2++; } //如果矩阵的秩小于增广矩阵的秩，则无解 if (rank1 < rank2) return null; //如果矩阵的秩小于方程组的数量，则有无穷解 if (rank1 < cols - 1) return new double[] { }; //从底往上依次求解 double[] result = new double[cols - 1]; for (int i = rows - 1; i >= 0; i--) { double y = matrix[i, cols - 1]; for (int j = i + 1; j < cols - 1; j++) { y -= matrix[i, j] * result[j]; } result[i] = matrix[i, i] != 0 ? Math.Round(y / matrix[i, i], 3) : 0; } return result; } /// <summary> /// 选主元。 /// </summary> /// <param name="matrix">增广系数矩阵。</param> /// <param name="i">当前行。</param> static void GaussPivoting(double[,] matrix, int i, int rows, int cols) { //选主元 int pivotRow = i; for (int k = i + 1; k < rows; k++) { if (Math.Abs(matrix[k, i]) > Math.Abs(matrix[pivotRow, i])) pivotRow = k; } //交换行 if (pivotRow != i) { double tmp = 0; for (int j = 0; j < cols; j++) { tmp = matrix[i, j]; matrix[i, j] = matrix[pivotRow, j]; matrix[pivotRow, j] = tmp; } } //除以主元，当前主元变为1 double pivot = matrix[i, i]; if (pivot != 0) { for (int j = i; j < cols; j++) matrix[i, j] /= pivot; } } }
【总结】
本文实现了两个算法，一是抛物线拟合，二是高斯消去法求解线性方程组。尤其是后者，是一个通用的算法，大家可以直接拿去使用。
不足之处是，本文没有考虑方程组无解情况下的抛物线处理，大家可以自己去完善一下，欢迎多交流。
【版权声明】
1、对于本文描述的算法，以及提供的源代码，保留所有权利。
2、本文的源代码只是供大家学习、研究之用，不得用于商业目的。
3、如果想在其它网站或平台转载，请征得本人同意。