C# 最小二乘法实现多点拟合圆心计算

在 C# 中，可以使用最小二乘法来实现多点拟合圆心计算。以下是一个示例代码：

using System;

namespace CircleFitting
{
    class Program
    {
        static void Main(string[] args)
        {
            // 输入多个点的坐标
            double[,] points = new double[,] 
            {
                {1, 1},
                {2, 3},
                {3, 5},
                {4, 7}
            };

            // 计算多点拟合圆心
            double centerX, centerY, radius;
            FitCircle(points, out centerX, out centerY, out radius);

            Console.WriteLine('圆心坐标：({0}, {1})', centerX, centerY);
            Console.WriteLine('半径：{0}', radius);

            Console.ReadLine();
        }

        static void FitCircle(double[,] points, out double centerX, out double centerY, out double radius)
        {
            int n = points.GetLength(0); // 点的数量

            // 构造系数矩阵A和常数向量b
            double[,] A = new double[n, 3];
            double[,] b = new double[n, 1];

            for (int i = 0; i < n; i++)
            {
                double x = points[i, 0];
                double y = points[i, 1];

                A[i, 0] = 2 * x;
                A[i, 1] = 2 * y;
                A[i, 2] = 1;

                b[i, 0] = x * x + y * y;
            }

            // 使用最小二乘法求解Ax=b的最小二乘解
            Matrix aMatrix = new Matrix(A);
            Matrix bMatrix = new Matrix(b);
            Matrix xMatrix = aMatrix.SolveLeastSquares(bMatrix);

            centerX = xMatrix[0, 0];
            centerY = xMatrix[1, 0];
            radius = Math.Sqrt(centerX * centerX + centerY * centerY + xMatrix[2, 0]);
        }
    }

    class Matrix
    {
        private double[,] data;
        public int RowCount { get; }
        public int ColumnCount { get; }

        public Matrix(double[,] data)
        {
            this.data = data;
            RowCount = data.GetLength(0);
            ColumnCount = data.GetLength(1);
        }

        public double this[int i, int j]
        {
            get { return data[i, j]; }
            set { data[i, j] = value; }
        }

        public Matrix SolveLeastSquares(Matrix b)
        {
            // 使用QR分解求解最小二乘解
            Matrix qr = QRDecomposition();
            Matrix q = qr.GetQMatrix();
            Matrix r = qr.GetRMatrix();

            Matrix bPrime = q.Transpose() * b;

            Matrix x = new Matrix(new double[ColumnCount, 1]);

            for (int i = ColumnCount - 1; i >= 0; i--)
            {
                double sum = 0;
                for (int j = i + 1; j < ColumnCount; j++)
                {
                    sum += r[i, j] * x[j, 0];
                }
                x[i, 0] = (bPrime[i, 0] - sum) / r[i, i];
            }

            return x;
        }

        private Matrix QRDecomposition()
        {
            Matrix q = new Matrix(new double[RowCount, ColumnCount]);
            Matrix r = new Matrix(new double[ColumnCount, ColumnCount]);

            Matrix a = Clone();
            for (int k = 0; k < ColumnCount; k++)
            {
                double norm = 0;
                for (int i = 0; i < RowCount; i++)
                {
                    norm += a[i, k] * a[i, k];
                }
                norm = Math.Sqrt(norm);

                r[k, k] = norm;

                for (int i = 0; i < RowCount; i++)
                {
                    q[i, k] = a[i, k] / norm;
                }

                for (int j = k + 1; j < ColumnCount; j++)
                {
                    double dotProduct = 0;
                    for (int i = 0; i < RowCount; i++)
                    {
                        dotProduct += q[i, k] * a[i, j];
                    }

                    r[k, j] = dotProduct;

                    for (int i = 0; i < RowCount; i++)
                    {
                        a[i, j] -= dotProduct * q[i, k];
                    }
                }
            }

            return new Matrix[] { q, r };
        }

        public Matrix Transpose()
        {
            double[,] transposedData = new double[ColumnCount, RowCount];
            for (int i = 0; i < ColumnCount; i++)
            {
                for (int j = 0; j < RowCount; j++)
                {
                    transposedData[i, j] = data[j, i];
                }
            }
            return new Matrix(transposedData);
        }

        public Matrix GetQMatrix()
        {
            Matrix qr = QRDecomposition();
            return qr[0];
        }

        public Matrix GetRMatrix()
        {
            Matrix qr = QRDecomposition();
            return qr[1];
        }

        public Matrix Clone()
        {
            double[,] clonedData = new double[RowCount, ColumnCount];
            Array.Copy(data, clonedData, data.Length);
            return new Matrix(clonedData);
        }
    }
}

在上面的代码中，我们首先定义了一个FitCircle方法，该方法接收一个多维数组points，其中存储了多个点的坐标。然后，我们使用最小二乘法来拟合圆心。

在FitCircle方法中，我们首先根据多个点的坐标构造了系数矩阵A和常数向量b。然后，我们使用Matrix类来实现最小二乘法，该类用于进行矩阵运算。我们使用QR分解来求解最小二乘解，然后计算出圆心坐标和半径。

最后，我们在Main方法中定义了一个示例点集，然后调用FitCircle方法来计算圆心坐标和半径，并将结果输出到控制台。

希望这个示例能帮助到你！