球面拟合算法

最新推荐文章于 2022-10-11 11:52:16 发布
最新推荐文章于 2022-10-11 11:52:16 发布
阅读量1.3w 收藏 112
点赞数 25
CC 4.0 BY-SA版权
分类专栏： 数值计算+数学软件
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/liyuanbhu/article/details/80201371

球面拟合算法

我在以前的博客中写过一篇介绍过最小二乘法圆拟合算法。前段时间有人给我留言，问我他的数据点是分布在一个球的球面上的，该如何拟合这个球的参数。球面拟合与圆的拟合其实是很类似的。只要按照我上篇文章的方法重新推导一下就能出来。正好这两天不太忙，在出差的路上推导了一下，把结果放到这里供大家参考。顺便多说一句，在高铁上用降噪耳机效果真的很好，很安静，值得拥有。

简单描述一下问题：我们有一组数据：(xi,yi,zi)(xi,yi,zi), 这些数据是分散在一个球面上的，当然数据中混杂着一些噪声。球面方程

(x−x0)2+(y−y0)2+(z−z0)2=R2(x−x0)2+(y−y0)2+(z−z0)2=R2

其中 (x0,y0,z0,R)(x0,y0,z0,R) 是我们需要求出的参数。

我们构造一个函数：

H(x0,y0,z0,R)=∑i=1N((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)2H(x0,y0,z0,R)=∑i=1N((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)2

我们求的 (x0,y0,z0,R)(x0,y0,z0,R) 是使 H 取得最小值的那组参数。也就是满足下面四个方程的(x0,y0,z0,R)(x0,y0,z0,R).
∂H∂x0=−4∑i=1N(xi−x0)((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∂H∂y0=−4∑i=1N(yi−z0)((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∂H∂z0=−4∑i=1N(zi−z0)((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∂H∂R=−4∑i=1NR((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∂H∂x0=−4∑i=1N(xi−x0)((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∂H∂y0=−4∑i=1N(yi−z0)((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∂H∂z0=−4∑i=1N(zi−z0)((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∂H∂R=−4∑i=1NR((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0

最后一个方程可以化简为：

∑i=1N((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∑i=1N((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0

利用这个式子可以把前三个方程化简为：
∑i=1Nxi((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∑i=1Nyi((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∑i=1Nzi((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∑i=1Nxi((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∑i=1Nyi((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0∑i=1Nzi((xi−x0)2+(yi−y0)2+(zi−z0)2−R2)=0

为了进一步化简，再做一次变换：
ui=xi−x¯vi=yi−y¯wi=zi−z¯ui=xi−x¯vi=yi−y¯wi=zi−z¯

这样替换之后有：
∑i=1Nui((ui−u0)2+(vi−v0)2+(wi−w0)2−R2)=0∑i=1Nvi((ui−u0)2+(vi−v0)2+(wi−w0)2−R2)=0∑i=1Nwi((ui−u0)2+(vi−v0)2+(wi−w0)2−R2)=0∑i=1N((ui−u0)2+(vi−v0)2+(wi−w0)2−R2)=0∑i=1Nui((ui−u0)2+(vi−v0)2+(wi−w0)2−R2)=0∑i=1Nvi((ui−u0)2+(vi−v0)2+(wi−w0)2−R2)=0∑i=1Nwi((ui−u0)2+(vi−v0)2+(wi−w0)2−R2)=0∑i=1N((ui−u0)2+(vi−v0)2+(wi−w0)2−R2)=0

再次化简后有：

(∑u2i)u0+(∑uivi)v0+(∑uiwi)w0=∑(u3i+uiv2i+uiw2i)2(∑uivi)u0+(∑v2i)v0+(∑viwi)w0=∑(u2ivi+v3i+viw2i)2(∑uiwi)u0+(∑viwi)v0+(∑w2i)w0=∑(u2iwi+v2iwi+w3i)2(∑ui2)u0+(∑uivi)v0+(∑uiwi)w0=∑(ui3+uivi2+uiwi2)2(∑uivi)u0+(∑vi2)v0+(∑viwi)w0=∑(ui2vi+vi3+viwi2)2(∑uiwi)u0+(∑viwi)v0+(∑wi2)w0=∑(ui2wi+vi2wi+wi3)2

解出这个方程就可以得到 (u0,v0,w0)(u0,v0,w0), 然后就可以得到 (x0,y0,z0)(x0,y0,z0) 了。之后带入第四个式子，就可以得到 RR 的值：
R=∑i=1N((ui−u0)2+(vi−v0)2+(wi−w0)2)N

下面是代码：

struct Point3D
{
    double x;
    double y;
    double z;
};

// 全选主元高斯消去法
//a-n*n 存放方程组的系数矩阵，返回时将被破坏
//b-常数向量
//x-返回方程组的解向量
//n-存放方程组的阶数
//返回0表示原方程组的系数矩阵奇异
int cagaus(double a[],double b[],int n,double x[])
{
    int *js, l, k, i, j, is, p, q;
    double d, t;
    js=new int [n];
    l=1;
    for (k=0;k<=n-2;k++)
    {
        d=0.0;
        for (i=k;i<=n-1;i++)
        {
            for (j=k;j<=n-1;j++)
            {
                t=fabs(a[i*n+j]);
                if (t>d)
                {
                    d=t;
                    js[k]=j;
                    is=i;
                }
            }
        }
        if (d+1.0==1.0)
        {
            l=0;
        }
        else
        {
            if (js[k]!=k)
            {
                for (i=0;i<=n-1;i++)
                {
                    p=i*n+k;
                    q=i*n+js[k];
                    t=a[p];
                    a[p]=a[q];
                    a[q]=t;
                }
            }
            if (is!=k)
            {
                for (j=k;j<=n-1;j++)
                {
                    p=k*n+j;
                    q=is*n+j;
                    t=a[p];
                    a[p]=a[q];
                    a[q]=t;
                }
                t=b[k];
                b[k]=b[is];
                b[is]=t;
            }
        }
        if (l==0)
        {
            delete js;
            return(0);
        }
        d=a[k*n+k];
        for (j=k+1;j<=n-1;j++)
        {
            p=k*n+j;
            a[p]=a[p]/d;
        }
        b[k]=b[k]/d;
        for (i=k+1;i<=n-1;i++)
        {
            for (j=k+1;j<=n-1;j++)
            {
                p=i*n+j;
                a[p]=a[p]-a[i*n+k]*a[k*n+j];
            }
            b[i]=b[i]-a[i*n+k]*b[k];
        }
    }
    d=a[(n-1)*n+n-1];
    if (fabs(d)+1.0==1.0)
    {
        delete js;
        return(0);
    }
    x[n-1]=b[n-1]/d;
    for (i=n-2;i>=0;i--)
    {
        t=0.0;
        for (j=i+1;j<=n-1;j++)
        {
            t=t+a[i*n+j]*x[j];
        }
        x[i]=b[i]-t;
    }
    js[n-1]=n-1;
    for (k=n-1;k>=0;k--)
    {
        if (js[k]!=k)
        {
            t=x[k];
            x[k]=x[js[k]];
            x[js[k]]=t;
        }
    }
    delete js;
    return 1;
}
bool sphereLeastFit(const std::vector<Point3D> &points,
                    double &center_x, double &center_y, double &center_z, double &radius)
{
    center_x = 0.0f;
    center_y = 0.0f;
    radius = 0.0f;

    int N = points.size();
    if (N < 4) // 最少 4 个点才能拟合出一个球面。
    {
        return false;
    }

    double sum_x = 0.0f, sum_y = 0.0f, sum_z = 0.0f;

    for (int i = 0; i < N; i++)
    {
        double x = points[i].x;
        double y = points[i].y;
        double z = points[i].z;
        sum_x += x;
        sum_y += y;
        sum_z += z;
    }
    double mean_x = sum_x / N;
    double mean_y = sum_y / N;
    double mean_z = sum_z / N;
    double sum_u2 = 0.0f, sum_v2 = 0.0f, sum_w2 = 0.0f;
    double sum_u3 = 0.0f, sum_v3 = 0.0f, sum_w3 = 0.0f;
    double sum_uv = 0.0f, sum_uw = 0.0f, sum_vw = 0.0f;
    double sum_u2v = 0.0f, sum_uv2 = 0.0f, sum_u2w = 0.0f, sum_v2w = 0.0f, sum_uw2 = 0.0f, sum_vw2 = 0.0f;
    for (int i = 0; i < N; i++)
    {
        double u = points[i].x - mean_x;
        double v = points[i].y - mean_y;
        double w = points[i].z - mean_z;
        double u2 = u * u;
        double v2 = v * v;
        double w2 = w * w;
        sum_u2 += u2;
        sum_v2 += v2;
        sum_w2 += w2;
        sum_u3 += u2 * u;
        sum_v3 += v2 * v;
        sum_w3 += w2 * w;
        sum_uv += u * v;
        sum_vw += v * w;
        sum_uw += u * w;
        sum_u2v += u2 * v;
        sum_u2w += u2 * w;
        sum_uv2 += u * v2;
        sum_v2w += v2 * w;
        sum_uw2 += u * w2;
        sum_vw2 += v * w2;
    }

    double A[3][3];
    double B[3] = {(sum_u3 + sum_uv2 + sum_uw2) / 2.0,
                   (sum_u2v + sum_v3 + sum_vw2) / 2.0,
                   (sum_u2w + sum_v2w + sum_w3) / 2.0};

    A[0][0] = sum_u2;
    A[0][1] = sum_uv;
    A[0][2] = sum_uw;
    A[1][0] = sum_uv;
    A[1][1] = sum_v2;
    A[1][2] = sum_vw;
    A[2][0] = sum_uw;
    A[2][1] = sum_vw;
    A[2][2] = sum_w2;

    double ans[3];

    if( cagaus((double *) A, B, 3, ans) == 0)
    {
        return false;
    }
    else
    {
        center_x = ans[0] + mean_x;
        center_y = ans[1] + mean_y;
        center_z = ans[2] + mean_z;

        double sum = 0;
        for(int i = 0; i < N; i++)
        {
            double dx = points[i].x - center_x;
            double dy = points[i].y - center_y;
            double dz = points[i].z - center_z;
            sum += dx * dx + dy * dy + dz * dz;
        }
        sum /= N;
        radius = sqrt(sum);
    }
    return true;
}
球面拟合-基于最小二乘法
huifer
04-11 5181
球面拟合-基于最小二乘法 球面公式 R2=x2+y2+z2 R^2 = x^2 + y^2+z^2 R2=x2+y2+z2 流程思路 最终结果是一个质地均匀的球体，并非椭圆，球面上的每一个点到圆心的距离相等 设球体中心为(x0,y0,z0)，球体表面点和中心点的关系如下 R2=(X−x0)2+(Y−y0)2+(Z−z0)2向量化简x2+y2+z2=2Xx0+2Yy0+2Zz0+r2−...
PCL：RANSAC球面拟合
孙悟空
12-16 3731
本文介绍了PCL进行RANSAC球面拟合的原理与方法。
球面数据拟合算法简介
热门推荐
Inspire
12-05 1万+
当我们手中握有大量的数据时，对于二维的数据，我们会对他们进行直线拟合、对数拟合，圆曲线的拟合等等。这些拟合的方法都是运用的了非常古老而又非常有效的方法，即最小二乘法。 今天给大家介绍一种三维球面数据的拟合方法，该方法也是运用的最小二乘的方法。旨在使拟合的半径在均方意义下误差达到最小。
c++ 球面拟合算法
08-03
球面拟合算法是一种用于确定一组离散点云数据在球面上的最佳拟合球的算法。该算法常被应用于三维重建、三维模型匹配、物体识别等领域。 球面拟合算法的基本思想是找到一个球面，使得该球面上的点到所有离散点的距离...
球面拟合算法 c++
最新发布
08-17
球面拟合算法是一种数学方法，用于找到一组数据点的最佳表面模型，通常假设这些数据点分布在球面上。在C++中，你可以通过优化技术，比如最小二乘法，来实现这个过程。一个常见的做法是使用多元回归或迭代算法，如...
非球面拟合球面的数据处理方法研究与实现.pdf
08-15
在现代光学加工领域，非球面拟合球面的数据处理方法扮演着至关重要的角色。随着科技的不断进步，对光学元件的精度要求越来越高，因此需要采用更加精细的检测方法和仪器来确保质量。轮廓仪作为精密测量仪器，能够提供...
利用多个二维点拟合圆、多个三维点拟合球面
12-21
总结来说，本主题涵盖了使用C++和最小二乘法进行数据拟合的基本概念，包括二维圆拟合和三维球面拟合的数学原理以及如何在实际编程中实现这些算法。这是一项实用的技能，尤其在处理几何数据、图像处理和物理建模等...
球面拟合
chj_kx_kx的博客
12-30 1464
使用Levenberg–Marquardt算法，参考https://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm。
RTK球面拟合
07-16
基于以上球心拟合的最小二乘构造原理，作者用MATLAB语言编写了球心拟合程序。为了估计拟合模型的合理性，作者利用数值模拟方式人为构造一组球数据，并用编写的拟合模型对球心进行拟合，拟合结果表明，模型最大拟合误差为优于1*10-4mm，因此球心的拟合模型是合理正确的。
空间球拟合程序（Matlab）
04-19
matlab程序 可以进行最小二乘法拟合空间中点坐标生成球方程
em算法matlab代码-Hypersphere-fitting:超球面拟合
05-26
em算法matlab代码超球面拟合 该存储库包含与我的超球面拟合工作相关的代码 IEEE标准 此文件夹包含与论文“ Julien Lesouple，Barbara Pilastre，Yoann Altmann和Jean-Yves Tourneret”已被IEEE信号处理快报所接受的论文“使用EM算法从噪声数据进行超球面拟合”相关的代码。 的MATLAB 此文件夹包含Matlab代码，该代码允许重现文章中的大多数结果，即图1和2。 图2、4、5和6以及表1，以及图1和2。 技术报告的2、3、4、5、6、7、9、10、11、12、13、14。 它借鉴了Matlab文件交换的两个功能： vmrand by Dylan Muir（）： Chen Yu-Hui（）的SphericalDistributionsRand： PYTHON 该文件夹包含一个轻巧的Jupyter笔记本，用于尝试提出的算法。
球体拟合（最小二乘法）：将球体拟合到一组噪声数据。 不需要宽弧或许多点。-matlab开发
05-30
给定一组数据点，该函数以最小二乘法计算数据的中心和半径。 最小二乘方程用于将矩阵求逆为 3x3，而不是直接在数据集上进行。 不需要大弧或许多数据点。 假设点不是奇异的（共面的）和真实的...... 在 R2010b 上创建，但应该适用于所有版本。
拟合圆算法
04-18
包括一个matlab代码，用来拟合程序图像中的圆，方便快捷。
最小二乘法球面拟合(附完整代码)
Marc Pony
12-04 7176
一种简洁、高效的球面拟合方法，附完整代码，方便学习
球面的最小二乘法拟合程序
weixin_41639107的博客
10-05 1317
%% 最小二乘法球面拟合程序 % 球面拟合算法 % H=load('X0.rpt');我在运行程序的时候用的是此行，但这样并不便于各位程序猿学习; % 为了便于各位程序猿学习，我将X0.RPT的内容列于下方，该数据文件来源于有限元分析软件PATRAN H= [97 39.657715 5.221072 0.199951 0.000000 -0...
高斯牛顿法进行三维球拟合(含C++代码)
jiangxing11的博客
10-11 2200
的最优解，最终求得极值，这里需要注意导数为0的点，不一定就是极值点，可能会是鞍点，也就是我们所说的局部最优解。然而，有些时候，若目标函数。太大，会导致其局部近似不精确，严重的时候，可能无法保证迭代收敛。下面介绍常见的两种优化算法：最速下降法(梯度下降法)，它是用于找到可微函数的局部最小值的一阶迭代优化算法，如果实值函数。当然，这样的算法也有一定的问题，例如：如果求出来的步长。当然，我们还需要该方向上取一个步长 λ，求得最快的下降方式，即。它的缺点是过于贪心，容易走出锯齿路线，反而增加了迭代次数。
梯度下降法拟合球面
sharemon的博客
06-17 744
今天遇到了球面拟合的问题，正好在上coursera上的机器学习的课程，干脆就用梯度下降拟合一个球面试试吧 %% clear clear; clc; close all; %% generate data theta = 0:2*pi/100:2*pi; beta = pi/6:2*pi/100:5*pi/6; r = 1; x = r*cos(beta)'*cos(theta) + 1; ...
评论 6
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值