ErrorAnalysis.git

			@@ -1,5 +1,6 @@
			using System;
			using System.Collections.Generic;
			using System.Drawing;
			using System.Linq;
			using System.Text;
			using System.Threading.Tasks;
			@@ -23,5 +24,322 @@
			double y = y0 + (value - x0) * (y1 - y0) / (x1 - x0);
			return y;
			}

			public static (double slope, double intercept) FitLine(PointF[] points)
			{
			// 计算各求和项
			int n = points.Length;
			double sumX = 0, sumY = 0, sumXY = 0, sumX2 = 0;

			foreach (var p in points)
			{
			sumX += p.X;
			sumY += p.Y;
			sumXY += p.X * p.Y;
			sumX2 += p.X * p.X;
			}

			// 计算斜率和截距
			double denominator = n * sumX2 - sumX * sumX;
			if (Math.Abs(denominator) < 1e-10) // 避免除以0
			throw new InvalidOperationException("点集过于集中，无法拟合直线");

			double slope = (n * sumXY - sumX * sumY) / denominator;
			double intercept = (sumY - slope * sumX) / n;

			return (slope, intercept);
			}

			public static double GetLineValue(double x, double slope, double intercept) => slope * x + intercept;

			public static (double a, double b, double c) FitParabola(PointF[] points)
			{
			if (points.Length < 3)
			throw new ArgumentException("至少需要3个点进行抛物线拟合");

			// 1. 计算各项求和值[1,6](@ref)
			int n = points.Length;
			double sx = 0, sy = 0, sx2 = 0, sx3 = 0, sx4 = 0, sxy = 0, sx2y = 0;

			foreach (var p in points)
			{
			double x = p.X;
			double y = p.Y;
			double x2 = x * x;
			double x3 = x2 * x;
			double x4 = x2 * x2;

			sx += x;
			sy += y;
			sx2 += x2;
			sx3 += x3;
			sx4 += x4;
			sxy += x * y;
			sx2y += x2 * y;
			}

			// 2. 构建正规方程矩阵[1,6](@ref)
			double[,] matrix = {
			{ sx4, sx3, sx2, sx2y },
			{ sx3, sx2, sx, sxy },
			{ sx2, sx, n, sy }
			};

			// 3. 使用高斯消元法求解
			double[] coefficients = GaussElimination(matrix);

			return (coefficients[0], coefficients[1], coefficients[2]);
			}

			// 高斯消元法实现（带部分主元选择）
			private static double[] GaussElimination(double[,] matrix)
			{
			int n = matrix.GetLength(0);
			double[] result = new double[n];

			// 前向消元
			for (int i = 0; i < n - 1; i++)
			{
			// 部分主元选择
			int maxRow = i;
			for (int k = i + 1; k < n; k++)
			{
			if (Math.Abs(matrix[k, i]) > Math.Abs(matrix[maxRow, i]))
			maxRow = k;
			}

			// 行交换
			if (maxRow != i)
			{
			for (int j = 0; j <= n; j++)
			{
			double temp = matrix[i, j];
			matrix[i, j] = matrix[maxRow, j];
			matrix[maxRow, j] = temp;
			}
			}

			// 消元过程
			for (int k = i + 1; k < n; k++)
			{
			double factor = matrix[k, i] / matrix[i, i];
			for (int j = i; j <= n; j++)
			{
			matrix[k, j] -= factor * matrix[i, j];
			}
			}
			}

			// 回代求解
			for (int i = n - 1; i >= 0; i--)
			{
			result[i] = matrix[i, n];
			for (int j = i + 1; j < n; j++)
			{
			result[i] -= matrix[i, j] * result[j];
			}
			result[i] /= matrix[i, i];
			}

			return result;
			}

			public static (double a, double b, double c) quadraticLine(double[] arrX, double[] arrY)
			{
			double[] coef = new double[3];
			double[,] AA = new double[3, 3];
			double[,] TAA = new double[3, 3];
			double[,] BAA = new double[3, 1];
			double[] B = new double[3];
			int m = arrX.Length;

			double[] arrXp = new double[m];
			double[] arrXc = new double[m];
			double[] arrXv = new double[m];
			double[] y = new double[m];

			for (int i = 0; i < m; i++)
			{
			arrXp[i] = arrX[i] * arrX[i];
			}
			for (int i = 0; i < m; i++)
			{
			arrXc[i] = arrX[i] * arrX[i] * arrX[i];
			}
			for (int i = 0; i < m; i++)
			{
			arrXv[i] = arrX[i] * arrX[i] * arrX[i] * arrX[i];
			}
			double a1 = mysum(arrX);
			double a2 = mysum(arrXp);
			double a3 = mysum(arrXc);
			double a4 = mysum(arrXv);

			AA[0, 0] = m;
			AA[0, 1] = a1;
			AA[0, 2] = a2;
			AA[1, 0] = a1;
			AA[1, 1] = a2;
			AA[1, 2] = a3;
			AA[2, 0] = a2;
			AA[2, 1] = a3;
			AA[2, 2] = a4;
			double[] xy = new double[m];
			double[] xxy = new double[m];
			for (int i = 0; i < m; i++)
			{
			xy[i] = arrX[i] * arrY[i];
			}
			for (int i = 0; i < m; i++)
			{
			xxy[i] = arrX[i] * arrX[i] * arrY[i];
			}
			double b1 = mysum(arrY);
			double b2 = mysum(xy);
			double b3 = mysum(xxy);
			B[0] = b1;
			B[1] = b2;
			B[2] = b3;
			TAA = ReverseMatrix(AA, 3);
			for (int i = 0; i < 3; i++)
			{
			for (int j = 0; j < 3; j++)
			{
			coef[i] = coef[i] + TAA[i, j] * B[j];
			}
			}
			return (coef[2], coef[1], coef[0]);

			}

			public static double[,] ReverseMatrix(double[,] dMatrix, int Level)
			{
			double dMatrixValue = MatrixValue(dMatrix, Level);
			if (dMatrixValue == 0) return null;
			double[,] dReverseMatrix = new double[Level, 2 * Level];
			double x, c;
			// Init Reverse matrix
			for (int i = 0; i < Level; i++)
			{
			for (int j = 0; j < 2 * Level; j++)
			{
			if (j < Level)
			dReverseMatrix[i, j] = dMatrix[i, j];
			else
			dReverseMatrix[i, j] = 0;
			}
			dReverseMatrix[i, Level + i] = 1;
			}
			for (int i = 0, j = 0; i < Level && j < Level; i++, j++)
			{
			if (dReverseMatrix[i, j] == 0)
			{
			int m = i;
			for (; dMatrix[m, j] == 0; m++) ;
			if (m == Level)
			return null;
			else
			{
			// Add i-row with m-row
			for (int n = j; n < 2 * Level; n++)
			dReverseMatrix[i, n] += dReverseMatrix[m, n];
			}
			}
			// Format the i-row with "1" start
			x = dReverseMatrix[i, j];
			if (x != 1)
			{
			for (int n = j; n < 2 * Level; n++)
			if (dReverseMatrix[i, n] != 0)
			dReverseMatrix[i, n] /= x;
			}
			// Set 0 to the current column in the rows after current row
			for (int s = Level - 1; s > i; s--)
			{
			x = dReverseMatrix[s, j];
			for (int t = j; t < 2 * Level; t++)
			dReverseMatrix[s, t] -= (dReverseMatrix[i, t] * x);
			}
			}
			// Format the first matrix into unit-matrix
			for (int i = Level - 2; i >= 0; i--)
			{
			for (int j = i + 1; j < Level; j++)
			if (dReverseMatrix[i, j] != 0)
			{
			c = dReverseMatrix[i, j];
			for (int n = j; n < 2 * Level; n++)
			dReverseMatrix[i, n] -= (c * dReverseMatrix[j, n]);
			}
			}
			double[,] dReturn = new double[Level, Level];
			for (int i = 0; i < Level; i++)
			for (int j = 0; j < Level; j++)
			dReturn[i, j] = dReverseMatrix[i, j + Level];
			return dReturn;
			}
			private static double MatrixValue(double[,] MatrixList, int Level)
			{
			double[,] dMatrix = new double[Level, Level];
			for (int i = 0; i < Level; i++)
			for (int j = 0; j < Level; j++)
			dMatrix[i, j] = MatrixList[i, j];
			double c, x;
			int k = 1;
			for (int i = 0, j = 0; i < Level && j < Level; i++, j++)
			{
			if (dMatrix[i, j] == 0)
			{
			int m = i;
			for (; dMatrix[m, j] == 0; m++) ;
			if (m == Level)
			return 0;
			else
			{
			// Row change between i-row and m-row
			for (int n = j; n < Level; n++)
			{
			c = dMatrix[i, n];
			dMatrix[i, n] = dMatrix[m, n];
			dMatrix[m, n] = c;
			}
			// Change value pre-value
			k *= (-1);
			}
			}
			// Set 0 to the current column in the rows after current row
			for (int s = Level - 1; s > i; s--)
			{
			x = dMatrix[s, j];
			for (int t = j; t < Level; t++)
			dMatrix[s, t] -= dMatrix[i, t] * (x / dMatrix[i, j]);
			}
			}
			double sn = 1;
			for (int i = 0; i < Level; i++)
			{
			if (dMatrix[i, i] != 0)
			sn *= dMatrix[i, i];
			else
			return 0;
			}
			return k * sn;
			}

			public static double mysum(double[] XX)
			{
			double reout = 0;
			int mm = XX.Length;
			for (int i = 0; i < mm; i++)
			{
			reout = reout + XX[i];
			}
			return reout;
			}

			public static double GetParabolaValue(double x, double a, double b, double c)
			{
			return a * x * x + b * x + c;
			}
			}
			}