using System;
|
using System.Collections.Generic;
|
using System.Drawing;
|
using System.Linq;
|
using System.Text;
|
using System.Threading.Tasks;
|
|
namespace ErrorAnalysis.Service
|
{
|
public class Utility
|
{
|
/// <summary>
|
/// 单线性插值法
|
/// </summary>
|
/// <param name="value">插值</param>
|
/// <param name="x0">起始值</param>
|
/// <param name="y0">起始值对应值</param>
|
/// <param name="x1">结束值</param>
|
/// <param name="y1">结束值对应值</param>
|
/// <returns>插值对应结果</returns>
|
public static double Interpolate(double value, double x0, double y0, double x1, double y1)
|
{
|
// Calculate the interpolated value using linear interpolation formula
|
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;
|
}
|
}
|
}
|
