c#最小二乘法拟合圆

Halcom

最小二乘法拟合圆
[1]https://blog.csdn.net/u014090257/article/details/108899982

[2]https://www.cnblogs.com/zxtang/p/15519549.html

1.         private void FittingCircle1D(List<System.Drawing.Point> pPointList, out double X_c, out double Y_c, out double R)
   
2.         {
   
3.             X_c = 0;
   
4.             Y_c = 0;
   
5.             R = 0;
   
6.             if (pPointList.Count < 3)
   
7.             {
   
8.                 return; // 最少3个点进行拟合
   
9.             }
   
10.             double X1 = 0;
    
11.             double Y1 = 0;
    
12.             double X2 = 0;
    
13.             double Y2 = 0;
    
14.             double X3 = 0;
    
15.             double Y3 = 0;
    
16.             double X1Y1 = 0;
    
17.             double X1Y2 = 0;
    
18.             double X2Y1 = 0;
    
19.             for (int i = 0; i < pPointList.Count; i++)
    
20.             {
    
21.                 X1 = X1 + pPointList[i].X;
    
22.                 Y1 = Y1 + pPointList[i].Y;
    
23.                 X2 = X2 + pPointList[i].X * pPointList[i].X;
    
24.                 Y2 = Y2 + pPointList[i].Y * pPointList[i].Y;
    
25.                 X3 = X3 + pPointList[i].X * pPointList[i].X * pPointList[i].X;
    
26.                 Y3 = Y3 + pPointList[i].Y * pPointList[i].Y * pPointList[i].Y;
    
27.                 X1Y1 = X1Y1 + pPointList[i].X * pPointList[i].Y;
    
28.                 X1Y2 = X1Y2 + pPointList[i].X * pPointList[i].Y * pPointList[i].Y;
    
29.                 X2Y1 = X2Y1 + pPointList[i].X * pPointList[i].X * pPointList[i].Y;
    
30.             }
    
31.             double C, D, E, G, H, N;
    
32.             double a, b, c;
    
33.             N = pPointList.Count;
    
34.             C = N * X2 - X1 * X1;
    
35.             D = N * X1Y1 - X1 * Y1;
    
36.             E = N * X3 + N * X1Y2 - (X2 + Y2) * X1;
    
37.             G = N * Y2 - Y1 * Y1;
    
38.             H = N * X2Y1 + N * Y3 - (X2 + Y2) * Y1;
    
39.             a = (H * D - E * G) / (C * G - D * D);
    
40.             b = (H * C - E * D) / (D * D - G * C);
    
41.             c = -(a * X1 + b * Y1 + X2 + Y2) / N;
    
42.             X_c = a / (-2);
    
43.             Y_c = b / (-2);
    
44.             R = Math.Sqrt(a * a + b * b - 4 * c) / 2;
    
45.         }
    
46.         

复制代码