【JS开源库】基于最小二乘法的离散点拟合圆形，计算圆心坐标和半径

介绍

使用的开源库：https://github.com/Meakk/circle-fit

一个用于快速拟合一组二维点的JavaScript库。
该实现基于最小二乘法来找到3个未知数：x和y坐标以及半径。

---

该算法不包含任何迭代求解器，将在线性时间内找到解，这意味着您将立即得到结果。
当您没有足够的非共线点通过返回错误标志来确定正确的圆时，此库会处理退化情况。

---

该库使用起来很简单。只有一个circlefit.js，全局对象circlefit会自动创建。此对象有3种方法：
addPoint(x,y)：在集合中添加一个点
clearPoints()：删除所有点
compute()：执行求解器并返回一个包含结果的对象。

---

返回的结果对象不仅包含圆的基本参数，还有其他参数：

• success（布尔值）：计算状态
• points（数组）：用户给出的所有点
• projections （数组）：将每个点投影到圆上的点
• distances （数组）：每个点到圆的距离
• center（Object）：圆的中心
• radius（数字）：圆的半径
• residue （数字）：最小二乘法的残差，可用于定义圆的质量
• computationTime（Number）：计算所花费的时间（以毫秒为单位）

示例

circlefit.js库的源码

出处：https://github.com/Meakk/circle-fit/blob/master/circlefit.js
是MIT协议。

源码复制到下面了：

var CIRCLEFIT = (function () {
  var my = {},
      points = [];

  function linearSolve2x2(matrix, vector) {
    var det = matrix[0]*matrix[3] - matrix[1]*matrix[2];
    if (det < 1e-8) return false; //no solution
    var y = (matrix[0]*vector[1] - matrix[2]*vector[0])/det;
    var x = (vector[0] - matrix[1]*y)/matrix[0];
    return [x,y];
  }

  my.addPoint = function (x, y) {
    points.push({x: x, y: y});
  }

  my.resetPoints = function () {
    points = [];
  }

  my.compute = function () {
    var result = {
      points: points,
      projections: [],
      distances: [],
      success: false,
      center: {x:0, y:0},
      radius: 0,
      residue: 0,
      computationTime: performance.now()
    };

    //means
    var m = points.reduce(function(p, c) {
      return {x: p.x + c.x/points.length,
              y: p.y + c.y/points.length};
    },{x:0, y:0});
    
    //centered points
    var u = points.map(function(e){
      return {x: e.x - m.x,
              y: e.y - m.y};
    });

    //solve linear equation
    var Sxx = u.reduce(function(p,c) {
      return p + c.x*c.x;
    },0);

    var Sxy = u.reduce(function(p,c) {
      return p + c.x*c.y;
    },0);

    var Syy = u.reduce(function(p,c) {
      return p + c.y*c.y;
    },0);

    var v1 = u.reduce(function(p,c) {
      return p + 0.5*(c.x*c.x*c.x + c.x*c.y*c.y);
    },0);

    var v2 = u.reduce(function(p,c) {
      return p + 0.5*(c.y*c.y*c.y + c.x*c.x*c.y);
    },0);

    var sol = linearSolve2x2([Sxx, Sxy, Sxy, Syy], [v1, v2]);

    if (sol === false) {
      //not enough points or points are colinears
      return result;
    }

    result.success = true;

    //compute radius from circle equation
    var radius2 = sol[0]*sol[0] + sol[1]*sol[1] + (Sxx+Syy)/points.length;
    result.radius = Math.sqrt(radius2);

    result.center.x = sol[0] + m.x;
    result.center.y = sol[1] + m.y;

    points.forEach(function(p) {
      var v = {x: p.x - result.center.x, y: p.y - result.center.y};
      var len2 = v.x*v.x + v.y*v.y;
      result.residue += radius2 - len2;
      var len = Math.sqrt(len2);
      result.distances.push(len - result.radius);
      result.projections.push({
        x: result.center.x + v.x*result.radius/len,
        y: result.center.y + v.y*result.radius/len
      });     
    });

    result.computationTime = performance.now() - result.computationTime;

    return result;
  }

  return my;
}());

//为了方便测试，笔者额外加的
module.exports = {
  CIRCLEFIT
}

测试用例1与效果图

var util=require("./circlefit.js");
var CIRCLEFIT = util.CIRCLEFIT;

var pts = [
    [137, 17672.52 ],
    [104.97, 17640.08],
    [12.52, 17539.34],
    [-45.99, 17469.43],
    [-101.92, 17397.44],
    [-180.79, 17285.75],
    [-253.44, 17169.91],
    [-319.65, 17050.28],
    [-379.22, 16927.2],
    [-431.98, 16801.05],
    [-463.28, 16715.43],
    [-491.44, 16628.72],
]

for (let i = 0; i < pts.length; i++) {
    CIRCLEFIT.addPoint(pts[i][0], pts[i][1]);
}

var result = CIRCLEFIT.compute();
if (result.success) {
    var c =  result.center;
    var x = Math.round(c.x);
    var y = Math.round(c.y);
    var r = Math.round(result.radius);
    console.log(`圆心=(${x},${y})`);
    console.log(`半径={${r}}`);
}

运行结果：

圆心=(1900,15900)
半径={2500}

效果图：

测试用例2与效果图

将点数组换成：

var pts = [
    [0, 0],
    [1, 1.1],
    [1.5, 0.8],
    [2, -0.5],
]

效果图：

补充

• 毕竟是最小二乘法来拟合，当点的噪音太大、太多时，可能拟合效果就不太好了。
• 可根据需要，转用其他语言重写。