最佳拟合线、圆和椭圆

• 下载演示可执行文件 - 2.6 MB
• 下载演示源代码 - 425.8 KB
• 仅下载库 - 15.5 KB

介绍 

该库用于确定一组输入点的最佳拟合二维线、圆或旋转椭圆。

以一组二维 x,y 点为例，这些点虽然近似圆，但并不精确。那么就存在一个中心点和半径，它们代表与这些点最匹配的圆。该库将为您提供该点和半径。类似地，对于一条线，您将获得斜率和 y 轴截距；对于一个椭圆，您将获得中心、长短半轴和旋转角度。

背景 

最佳拟合形状在工程、测量和计算几何学中有应用。这方面的文献很多，但可重用的代码似乎很少。该库提供了一种将此功能集成到客户端应用程序中的方法。它使用最小二乘法来实现这一点。

使用代码

代码是标准的 C++，依赖于 Boost ublas 库进行矩阵运算。提供了 Visual Studio 2008 解决方案文件和 Makefile，分别用于在 Windows 和 Linux 上构建。Boost 需要预先安装，并在您的构建环境中设置相应的目录路径。

总体概述

• 使用工厂类创建 BestFit 派生计算对象的实例。
• 决定要将计算的文本输出发送到哪个 std::ostream（如果有）。例如，输出到控制台或 std::stringstream，或者完全不输出。
• 创建 BestFitIO 结构的两个实例，一个用于输入参数，一个用于输出。
• 至少需要设置输入结构中的两个成员：numPoints 是要提供的点数，points 是一个双精度浮点数对序列的起始地址，表示 x,y 点。
• 输出结构可以保留为默认构造，但请参阅下面的 进一步的输入和输出 部分以获取其他选项。
• 在计算对象上调用 Compute()，传入 I/O 结构。
• 如果需要，从输出 BestFitIO 实例获取结果，或仅依赖于输出到 ostream 的内容。
• delete 工厂创建的对象。

一个简单的例子：

这是一个计算 10 个测试点的最佳拟合圆的示例，这些点近似半径为 1、中心为 (100,200) 的圆。然后发现这些点的真实中心为 (99.969,200.016)，半径为 0.977。

#include <iostream>
#include <sstream>
#include "BestFit.h"

int main(int argc, char * argv[])
{
    int count = 10;
    double points[20] = { 
        99.9510022,201.01559,
        100.0701559,199.016258,
        100.921381,200.0685969,
        99.031849,199.895546,
        99.469357,200.869111,
        100.525856,200.793908,
        99.13055,200.56031,
        100.787632,199.483885,
        100.498142,199.169206,
        99.363449,199.276086
    };

    // Input settings specify the source coordinates.
    BestFitIO input;
    input.numPoints = count;
    input.points = points;

    // Output settings in this simple example left with default values.
    BestFitIO output;

    // Create the best-fitting circle object, and make it do the work.
    int type = BestFitFactory::Circle;
    BestFit *b = BestFitFactory::Create(type, std::cout);
    b->Compute(input, output);

    ... // output contains solved parameters, and optionally the adjusted points.

    // The factory allocated memory needs to be freed.  
    delete b;
    return 0;
}

获取结果

计算后，BestFitIO::outputFields 包含形状的参数，所有参数都是 double 类型。以下枚举用于索引此数组

enum { LineGradient, LineYIntercept };
enum { CircleCentreX, CircleCentreY, CircleRadius };
enum { EllipseCentreX, EllipseCentreY, EllipseMajor, EllipseMinor, EllipseRotation };

例如，要检索上面示例中圆的已解算参数

double centre_x = output.outputFields[BestFitIO::CircleCentreX];
double centre_y = output.outputFields[BestFitIO::CircleCentreY];
double radius   = output.outputFields[BestFitIO::CircleRadius];

此外，一些质量统计信息会输出到指定的 ostream，在本例中是 std::cout。通过输入 BestFitIO::verbosity 设置可以设置输出量的大小。通过工厂构造函数的单参数重载也可以实现无输出。

Global check of the adjustment     ***PASSES***
Check of the evaluated unknowns    ***PASSES***
Number of observations             10
Number of unknowns                 3
Degrees of freedom                 7
Iterations until convergence       6
Variance                           0.000690
Std. dev. observation unit weight  0.026277
***********************************
Adj. circle centre X               99.968823
Adj. circle centre Y               200.016389
Adj. circle radius                 0.977187

这是一个圆作为最佳拟合形状的示例 - 对于线或椭圆，遵循相同的过程。

进一步的输入和输出

除了求解几何参数外，该库还可以计算

1. 调整后的点 - 对于每个输入点，对应于最佳拟合形状的点
2. 残差 - 对输入点的修正

在上面的简单示例中，没有输出调整后的点。但是，如果 output.wantAdjustedObs=true，则 output.points 数组将包含上面 (1) 中描述的点。可以为 output.points 分配内存，或者为了节省内存，可以将 output.points = input.points。然后输入点将被调整后的点覆盖。另外，通过设置 output.wantResiduals=true，将分配 output.residuals 数组并填充上面 (2) 中描述的修正。这两个谓词的默认值都为 false。 

BestFitIO output;
output.points = new double[20]; // allocate space for adjusted points.
output.wantAdjustedObs = true;  // output.points will get populated
output.wantResiduals = true;    // output.residuals will get allocated and populated

if (bestFit.Compute(input, output))
{
    ... // extract data from output.points and output.residuals
    delete [] output.residuals;
    delete [] output.points;
}

演示应用程序

作为一个测试工具和演示该库使用的图形化方法，我编写了一个演示应用程序。它依赖于 MFC 和 Dundas Ultimate Grid，因此仅适用于 Windows。请注意，通过使用该库，它也依赖于 Boost。

演示应用程序将最佳拟合代码用作静态链接库，并允许您

• 查看原始和调整后的坐标数据，以及点的图形视图；
• 右键单击添加点，左键单击从图形视图中删除点；
• 动态查看调整结果的变化；
• 打开和保存文件中的点；
• 生成近似给定线、圆或椭圆参数的随机分布点的测试数据。

关注点

命令行程序

库代码也可以编译为独立的控制台程序。它从 stdio 接收流输入。要这样构建，您只需要在预处理器定义中取消定义 BEST_FIT_LIB 并将构建目标设置为 .exe。对于 Linux，Makefile 默认目标就是此模式。当定义 BEST_FIT_LIB 构建时，有一个额外的依赖项 Boost program_options 用于解析命令行开关。命令行用法如下

$ ./best-fit
Usage: best-fit [options] NUM_POINTS X,Y ...
Allowed options:
  -t [ --type ] arg     type of equation (0=line, 1=circle, 2=ellipse)
  -v [ --verbose ] arg  0=simple, 1=default, 2=verbose, 3=more verbose, etc.
  -h [ --help ]         produce this message

速度/性能指标 

在 Core 2 Duo 上测试最佳拟合的速度，对于增加的点数，显示以下速度统计信息。

	未知数数量	100 个点	1000 个点	10000 个点
Line	2	<1毫秒	2毫秒	15毫秒
圆	3	<1毫秒	4毫秒	60毫秒
椭圆	5	1毫秒	10毫秒	124毫秒

Release 和 Debug 版本的 Boost uBlas 在矩阵乘法时间上存在巨大的差异。这在互联网上已有充分记载。我在本文章的早期版本中不幸忽略了这一点，该版本显示了严重不正确的时序。这已得到纠正，上述时序是在定义 NDEBUG 以禁用耗时的调试模式检查的情况下进行的。

算法 

线性最小二乘调整是代码背后的算法。具体来说，是调整的参数化情况。给定未知数的初始值，这种最小二乘法会求解小的调整，这些调整会优化初始值，然后迭代使用，直到解决方案收敛。矩阵运算是其中不可或缺的一部分，不幸的是包括矩阵求逆，这意味着对于大型矩阵它会很慢。

在最佳拟合线的情况下，只有 y 值被视为观测值。x 值被视为不变。最小化每个点到直线的垂直最短距离的平方和。这是最小二乘调整的普通参数化情况。

在圆和旋转椭圆的情况下，x 和 y 值都被视为观测值。最小化每个点到形状的距离的平方和。这是通用最小二乘调整的一个特例，称为拟参数化情况。

参数化和拟参数化情况都以相同的方式求解。区别在于权重矩阵的构造以及残差的含义。拟参数化情况的残差不像参数化情况那样直接等同于有用的几何量。

解向量（或未知数的修正）由下式给出

残差向量（或观测值的修正）由下式给出

其中 A 是观测值的设计矩阵，P 是权重矩阵，ℓ 是观测向量。解向量方程中存在的逆运算需要矩阵求逆，并对性能产生负面影响。

为了形成设计矩阵，几何形状的方程需要是线性的。直线方程已经是线性的，但圆和旋转椭圆的方程需要先通过泰勒展开进行线性化。

首先通过近似确定未知数的初始值。将每个坐标观测值和初始未知值代入线性化方程，形成 A 矩阵的一行。ℓ 矩阵是观测值减去初始值。

然后确定解向量 x，并将其添加到初始值中以进行优化。如果变化足够小，则解决方案已收敛，初始值就是最终值。否则，新的初始值将用于重新构建矩阵 A、P 和 ℓ，然后再进行一次迭代。

总之，这是一种方法，不是唯一的方法，也不是最好的方法。然而，这种方法不需要考虑直线、圆和椭圆之间的差异，这使得它非常方便，可以简洁地实现为一个通用库。

致谢和来源：Heinz Rüther，高级工程测量讲义（未发表），地理空间信息系，开普敦大学，1996 年。

历史

• 2013 年 12 月：首次发布；
• 2014 年 1 月：根据评论中的反馈添加了算法部分。
• 2014 年 2 月：修复了几个错误 修复 和改进了测试套件 改进，将 Makefile 默认更改为 Release 构建，但最重要的是更新了时序部分。
• 2017 年 1 月：根据评论中的反馈，添加了演示应用程序作为预编译的可执行文件，方便下载。