zoukankan      html  css  js  c++  java
  • 最小二乘法拟合圆

     

     

     

     好了。下面给出个代码,这个代码的具体公式和我这里给出的有一点小差异,但是原理是相同的。

     1 /**
     2  * 最小二乘法拟合圆
     3  * 拟合出的圆以圆心坐标和半径的形式表示
     4  * 此代码改编自 newsmth.net 的 jingxing 在 Graphics 版贴出的代码。
     5  * 版权归 jingxing, 我只是搬运工外加一些简单的修改工作。
     6  */
     7 typedef complex<int> POINT;
     8 bool circleLeastFit(const std::vector<POINT> &points, double &center_x, double &center_y, double &radius)
     9 {
    10      center_x = 0.0f;
    11      center_y = 0.0f;
    12      radius = 0.0f;
    13      if (points.size() < 3)
    14      {
    15          return false;
    16      }
    17 
    18      double sum_x = 0.0f, sum_y = 0.0f;
    19      double sum_x2 = 0.0f, sum_y2 = 0.0f;
    20      double sum_x3 = 0.0f, sum_y3 = 0.0f;
    21      double sum_xy = 0.0f, sum_x1y2 = 0.0f, sum_x2y1 = 0.0f;
    22 
    23      int N = points.size();
    24      for (int i = 0; i < N; i++)
    25      {
    26          double x = points[i].real();
    27          double y = points[i].imag();
    28          double x2 = x * x;
    29          double y2 = y * y;
    30          sum_x += x;
    31          sum_y += y;
    32          sum_x2 += x2;
    33          sum_y2 += y2;
    34          sum_x3 += x2 * x;
    35          sum_y3 += y2 * y;
    36          sum_xy += x * y;
    37          sum_x1y2 += x * y2;
    38          sum_x2y1 += x2 * y;
    39      }
    40 
    41      double C, D, E, G, H;
    42      double a, b, c;
    43 
    44      C = N * sum_x2 - sum_x * sum_x;
    45      D = N * sum_xy - sum_x * sum_y;
    46      E = N * sum_x3 + N * sum_x1y2 - (sum_x2 + sum_y2) * sum_x;
    47      G = N * sum_y2 - sum_y * sum_y;
    48      H = N * sum_x2y1 + N * sum_y3 - (sum_x2 + sum_y2) * sum_y;
    49      a = (H * D - E * G) / (C * G - D * D);
    50      b = (H * C - E * D) / (D * D - G * C);
    51      c = -(a * sum_x + b * sum_y + sum_x2 + sum_y2) / N;
    52 
    53      center_x = a / (-2);
    54      center_y = b / (-2);
    55      radius = sqrt(a * a + b * b - 4 * c) / 2;
    56      return true;
    57 }

    下图是个实际测试的结果。对这种均匀分布的噪声数据计算的结果还是很准确的。

    这里写图片描述

    但是当数据中有部分偏向某一个方向的干扰数据时。结果就不那么乐观了。下图就很说明问题。

    这里写图片描述

    数据点中有 20% 是干扰数据。拟合出的圆就明显被拽偏了。

     

  • 相关阅读:
    XAML学习笔记之Layout(五)——ViewBox
    XAML学习笔记——Layout(三)
    XAML学习笔记——Layout(二)
    XAML学习笔记——Layout(一)
    从0开始搭建SQL Server 2012 AlwaysOn 第三篇(安装数据,配置AlwaysOn)
    从0开始搭建SQL Server 2012 AlwaysOn 第二篇(配置故障转移集群)
    从0开始搭建SQL Server 2012 AlwaysOn 第一篇(AD域与DNS)
    Sql Server 2012 事务复制遇到的问题及解决方式
    Sql Server 2008R2升级 Sql Server 2012 问题
    第一次ACM
  • 原文地址:https://www.cnblogs.com/ybqjymy/p/14201802.html
Copyright © 2011-2022 走看看