切比雪夫(最小区域法)平面拟合算法

这篇具有很好参考价值的文章主要介绍了切比雪夫(最小区域法)平面拟合算法。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

欢迎关注更多精彩
关注我,学习常用算法与数据结构,一题多解,降维打击。

本期话题:切比雪夫(最小区域法)平面拟合算法

相关背景和理论
点击前往
主要介绍了应用背景和如何转化成线性规划问题

切比雪夫拟合,高阶算法,c++,数学,切比雪夫拟合,最小区域法,平面拟合,线性规划

平拟合输入和输出要求

输入

  1. 10到631个点,全部采样自平面附近。
  2. 每个点3个坐标,坐标精确到小数点后面20位。
  3. 坐标单位是mm, 范围[-500mm, 500mm]。

输出

  1. 平面上一点X0,用三个坐标表示。
  2. 法向A。
  3. 平面度F,所有点到平面距离最大的2倍。

精度要求

  1. X0点平面距离不能超过0.0001mm。
  2. 法向A与标准法向夹角不能超过0.0000001rad。
  3. F与标准平面度误差不能超过0.00001mm。

平面优化标函数

根据认证要求,平面拟合转化成数学表示如下:

平面参数化表示

  1. 平面上1点X0 = (x0, y0, z0)。
  2. 方向单位向量A=(a,b,c)。

点到平面距离

第i个点 pi(xi, yi, zi)。

根据点乘得到距离

d i = ∥ ( p i − X 0 ) ⋅ A ∥ ∥ A ∥ d_i = \frac { \left \| (p_i-X_0)\cdot A \right \|}{\left \| A \right \|} di=A(piX0)A

展开一下:

d i = a ( x i − x 0 ) + b ( y i − y 0 ) + c ( z i − z 0 ) a 2 + b 2 + c 2 d_i = \frac {a(x_i-x_0) + b(y_i-y_0) + c(z_i-z_0)} {\sqrt{a^2+b^2+c^2}} di=a2+b2+c2 a(xix0)+b(yiy0)+c(ziz0)

优化能量方程

能量方程 H = f ( X 0 , A ) = max ⁡ 1 n ∣ d i ∣ H=f(X0, A)=\displaystyle \max_1^n {|d_i|} H=f(X0,A)=1maxndi

X0, A是未知量,拟合平面的过程也可以理解为优化X0, A使得方程H最小。

如何3个参数表示平面

如果直接拿6个参数去做迭代,1是比较麻烦,会出现比较难解的方向,2是平面上的点有很多个,结果不唯一。

当平面法向与Z轴偏差比较小的时候可以使用3个参数来表示。

切比雪夫拟合,高阶算法,c++,数学,切比雪夫拟合,最小区域法,平面拟合,线性规划

如上图,绿线为Z轴,橙色线为XOY平面。

由于法向与Z轴比较相近,可以设法向为(a, b, 1), a,b 是比较小的量。

现在表示 平面还需要一个点,规定点必须在以(a,b,1)为法向过0点的直线上。

就有直线公式 x0/a=y0/b=z0/c

那么只要知道z0就可知 x0=az0, y0=bz0.

综上,可以使用a, b, z0 表示一个法向与Z轴相近的 平面。

转化为线性规划(法向与Z轴接近)

设 a = ( z 0 , A a , A b ) , d i = F ( x i ;   a ) , 引入 Γ = M A X i = 1 n    ∣ d i ∣ 设a=(z_0, A_a, A_b), d_i=F(x_i;\ a), 引入\Gamma=\overset n{\underset {i=1}{MAX}}\;|d_i| a=(z0,Aa,Ab),di=F(xi; a),引入Γi=1MAXndi

根据上述定义,可以将原来的最值问题转化为下述条件

对于所有点应该满足

F ( x i ;   a ) ≤ Γ , ( F ( x i ;   a ) > 0 ) F(x_i;\ a)\le \Gamma, (F(x_i;\ a)>0) F(xi; a)Γ,(F(xi; a)>0)

− F ( x i ;   a ) ≤ Γ , ( F ( x i ;   a ) < 0 ) -F(x_i;\ a)\le \Gamma, (F(x_i;\ a)<0) F(xi; a)Γ,(F(xi; a)<0)

我们可以通过小量迭代慢慢减小Γ

m a x      Δ Γ s . t .     F ( x i , a ) + J Δ a ≤ Γ − Δ Γ , ( i = 1 , 2... n )           − ( F ( x i , a ) + J Δ a ) ≤ Γ − Δ Γ , ( i = 1 , 2... n ) Δ Γ ≥ 0 \begin {array}{c}max \ \ \ \ \Delta {\Gamma}\\ s.t.\ \ \ F(x_i, a) + J\Delta a \le \Gamma -\Delta \Gamma, (i=1,2...n)\\ \ \ \ \ \ \ \ \ \ -(F(x_i, a) + J\Delta a) \le \Gamma -\Delta \Gamma, (i=1,2...n)\\ \Delta \Gamma \ge0\end{array} max    ΔΓs.t.   F(xi,a)+JΔaΓΔΓ,(i=1,2...n)         (F(xi,a)+JΔa)ΓΔΓ,(i=1,2...n)ΔΓ0

上述条件不需要管 F ( x i , a ) + J Δ a 正负情况,若 F ( x i , a ) + J Δ a 为正 − ( F ( x i , a ) + J Δ a ) ≤ Γ − Δ Γ 必成立,反之亦然。 上述条件不需要管F(x_i, a) + J\Delta a正负情况,若F(x_i, a) + J\Delta a为正-(F(x_i, a) + J\Delta a) \le \Gamma -\Delta \Gamma必成立,反之亦然。 上述条件不需要管F(xi,a)+JΔa正负情况,若F(xi,a)+JΔa为正(F(xi,a)+JΔa)ΓΔΓ必成立,反之亦然。
求解出以后更新a, Γ。

对线性规划模型标准化

需要对 Δ z 0 , Δ A a , Δ A b 拆解,要求变量都要大于等于 0 需要对\Delta z_0, \Delta A_a, \Delta A_b 拆解,要求变量都要大于等于0 需要对Δz0,ΔAa,ΔAb拆解,要求变量都要大于等于0

m a x      Δ Γ s . t .     J i ⋅ [ Δ z 0 + - Δ z 0 - Δ A a + - Δ A a - Δ A b + − Δ A b − ] + Δ Γ ≤ Γ - d i , ( i = 1 , 2... n )        − J i ⋅ [ Δ z 0 + - Δ z 0 - Δ A a + - Δ A a - Δ A b + − Δ A b − ] + Δ Γ ≤ Γ + d i , ( i = 1 , 2... n ) Δ z 0 + , Δ z 0 − , Δ A a + , Δ A a − , Δ A b + , Δ A b − , Δ Γ ≥ 0 ( 1 ) \begin {array}{c}max \ \ \ \ \Delta {\Gamma}\\ s.t.\ \ \ J_i \cdot \begin {bmatrix} \Delta z_0^+-\Delta z_0^-\\ \Delta A_a^+-\Delta A_a^-\\ \Delta A_b^+- \Delta A_b^- \end{bmatrix} +\Delta \Gamma\le \Gamma-d_i, (i=1,2...n)\\\\ \ \ \ \ \ \ -J_i \cdot \begin {bmatrix} \Delta z_0^+-\Delta z_0^-\\ \Delta A_a^+-\Delta A_a^-\\ \Delta A_b^+- \Delta A_b^- \end{bmatrix}+\Delta \Gamma\le \Gamma+d_i, (i=1,2...n)\\ \Delta z_0^+, \Delta z_0^-, \Delta A_a^+, \Delta A_a^-, \Delta A_b^+, \Delta A_b^-, \Delta \Gamma \ge0\end{array} (1) max    ΔΓs.t.   Ji Δz0+Δz0ΔAa+ΔAaΔAb+ΔAb +ΔΓΓdi,(i=1,2...n)      Ji Δz0+Δz0ΔAa+ΔAaΔAb+ΔAb +ΔΓΓ+di,(i=1,2...n)Δz0+,Δz0,ΔAa+,ΔAa,ΔAb+,ΔAb,ΔΓ0(1)

算法描述

法向与Z轴重合时
x 0 = 0 , y 0 = 0 , z 0 = 0 , a = 0 , b = 0 , c = 1 x_0=0, y_0=0, z_0=0, a=0,b =0, c=1 x0=0,y0=0,z0=0,a=0,b=0,c=1

d i = a ( x i − x 0 ) + b ( y i − y 0 ) + c ( z i − z 0 ) a 2 + b 2 + c 2 = z i d_i = \frac {a(x_i-x_0) + b(y_i-y_0) + c(z_i-z_0)} {\sqrt{a^2+b^2+c^2}}=z_i di=a2+b2+c2 a(xix0)+b(yiy0)+c(ziz0)=zi

J, D的计算。

3个未知分别对d_i求导结果如下:

求导后a=b=z0=x0=y0=0,c=1要代入

∂ d i ∂ z 0 = ∂ a ( x i − x 0 ) + b ( y i − y 0 ) + c ( z i − z 0 ) a 2 + b 2 + c 2 ∂ z 0 ∣ z 0 = 0 = − 1 \frac {\partial d_i} {\partial z0}=\frac {\partial \frac {a(x_i-x_0) + b(y_i-y_0) + c(z_i-z_0)}{\sqrt{a^2+b^2+c^2}}} {\partial z0}_{|z0=0} = -1 z0di=z0a2+b2+c2 a(xix0)+b(yiy0)+c(ziz0)z0=0=1

∂ d i ∂ a = ∂ a ( x i − x 0 ) + b ( y i − y 0 ) + c ( z i − z 0 ) a 2 + b 2 + c 2 ∂ a ∣ a = 0 = x i \frac {\partial d_i} {\partial a}=\frac {\partial \frac {a(x_i-x_0) + b(y_i-y_0) + c(z_i-z_0)}{\sqrt{a^2+b^2+c^2}}} {\partial a}_{|a=0} = x_i adi=aa2+b2+c2 a(xix0)+b(yiy0)+c(ziz0)a=0=xi

∂ d i ∂ b = ∂ a ( x i − x 0 ) + b ( y i − y 0 ) + c ( z i − z 0 ) a 2 + b 2 + c 2 ∂ b ∣ b = 0 = y i \frac {\partial d_i} {\partial b}=\frac {\partial \frac {a(x_i-x_0) + b(y_i-y_0) + c(z_i-z_0)}{\sqrt{a^2+b^2+c^2}}} {\partial b}_{|b=0} = y_i bdi=ba2+b2+c2 a(xix0)+b(yiy0)+c(ziz0)b=0=yi

一次迭代过程

  1. 确定平面初值,Γ, 任意选择三点拟合平面

  2. 根据上述公式(1)构建线性规划方程

  3. 求解 Δ p \Delta p Δp

  4. 更新解
    [ x 0 y 0 z 0 ] = [ x 0 y 0 z 0 ] + U T ⋅ [ p a ∗ p z 0 p b ∗ p z 0 p z 0 ] [ a b c ] = U T ⋅ [ p a p b 1 ] . n o r m a l i z e ( ) Γ = Γ − Δ Γ \begin {array}{l} \begin {bmatrix}x_0 \\ y_0 \\ z_0 \end {bmatrix} = \begin {bmatrix}x_0 \\ y_0 \\ z_0 \end {bmatrix} + U^T \cdot \begin{bmatrix}p_a*p_{z_0} \\ p_b*p_{z_0} \\ p_{z_0} \end {bmatrix}\\ \\ \begin {bmatrix}a \\ b \\ c \end {bmatrix} = U^T \cdot \begin{bmatrix}p_a \\ p_b \\ 1 \end {bmatrix}.normalize() \\ \\ \Gamma = \Gamma-\Delta \Gamma \end {array} x0y0z0 = x0y0z0 +UT papz0pbpz0pz0 abc =UT papb1 .normalize()Γ=ΓΔΓ

  5. 重复2直到收敛

最后,输出时F=2*Γ

代码实现

代码链接:https://gitcode.com/chenbb1989/3DAlgorithm/blob/master/CBB3DAlgorithm/Fitting/chebyshev/PlaneFitter.cpp

拟合代码

#include "PlaneFitter.h"
#include "../gauss/FittingPlane.h"
#include <algorithm>
#include <Eigen/Dense>


namespace Chebyshev {
	double PlaneFitter::F(Fitting::Plane plane, const Point& p)
	{
		auto de = p - plane.BasePoint;
		return plane.Orient.dot(de);
	}

	double PlaneFitter::GetError(Fitting::Plane plane, const std::vector<Eigen::Vector3d>& points)
	{
		double err = 0;
		for (auto& p : points) { 
			err = std::max(err, abs(F(plane, p)));
		}

		return err;
	}

	Fitting::Matrix PlaneFitter::Jacobi(const std::vector<Eigen::Vector3d>& points)
	{
		Fitting::Matrix J(points.size(), 3);
		for (int i = 0; i < points.size(); ++i) {
			auto& p = points[i];
			J(i, 0) = -1;
			J(i, 1) = p.x();
			J(i, 2) = p.y();
		}
		return J;
	}

	void PlaneFitter::beforHook(const std::vector<Eigen::Vector3d>& points)
	{
		U = Fitting::getRotationByOrient(plane.Orient);
		for (int i = 0; i < points.size(); ++i) {
			transPoints[i] = U * (points[i] - plane.BasePoint);
		}
	}

	void PlaneFitter::afterHook(const Eigen::VectorXd& xp)
	{
		plane.BasePoint += U.transpose() * Eigen::Vector3d(xp(0) * xp(1), xp(0) * xp(2), xp(0));
		plane.Orient = U.transpose() * Eigen::Vector3d(xp(1), xp(2), 1).normalized();
		gamma -= xp(3);
	}
	Eigen::VectorXd PlaneFitter::getDArray(const std::vector<Eigen::Vector3d>& points)
	{
		Eigen::VectorXd D(points.size());
		for (int i = 0; i < points.size(); ++i)D(i) = points[i].z();
		return D;
	}
	bool PlaneFitter::GetInitFit(const std::vector<Eigen::Vector3d>& points)
	{
		if (points.size() < 3)return false;
		Fitting::FittingBase* fb = new Gauss::FittingPlane();
		fb->Fitting(points, &plane);
		delete fb;

		gamma = GetError(plane, points);
		return true;
	}
	double PlaneFitter::F(const Eigen::Vector3d& p)
	{
		return Chebyshev::PlaneFitter::F(plane, p);
	}
	double PlaneFitter::GetError(const std::vector<Eigen::Vector3d>& points)
	{
		return Chebyshev::PlaneFitter::GetError(plane, points);
	}
	void PlaneFitter::Copy(void* ele)
	{
		memcpy(ele, &plane, sizeof(Fitting::Plane));
	}
	PlaneFitter::PlaneFitter()
	{
		ft = Fitting::FittingType::CHEBYSHEV;
	}
}


测试结果

https://gitcode.com/chenbb1989/3DAlgorithm/blob/master/CBB3DAlgorithm/Fitting/chebyshev/chebyshev-testdata/officialtest/fitting_result/result.txt
C17 : PLANE : pass
C18 : PLANE : pass
C19 : PLANE : pass
C20 : PLANE : pass
C21 : PLANE : pass
C22 : PLANE : pass
C23 : PLANE : pass
C24 : PLANE : pass
C25 : PLANE : pass
C26 : PLANE : pass


本人码农,希望通过自己的分享,让大家更容易学懂计算机知识。创作不易,帮忙点击公众号的链接。

切比雪夫拟合,高阶算法,c++,数学,切比雪夫拟合,最小区域法,平面拟合,线性规划文章来源地址https://www.toymoban.com/news/detail-843958.html

到了这里,关于切比雪夫(最小区域法)平面拟合算法的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 切比雪夫不等式,大数定律及极限定理。

    1.定理 若随机变量X的期望EX和方差DX存在,则对任意ε 0,有    P{ |X - EX| = ε } = DX/ε 2 或 P{ |X - EX| ε } = 1 - DX/ε 2 2.解析定理 ①该定理对 X 服从什么分布不做要求,仅EX DX存在即可。 ②“| |” 由于X某次试验结果可能大于期望值,也可能小于期望值,但总在其旁边波动,所 以加

    2024年02月06日
    浏览(36)
  • 马尔可夫不等式、切比雪夫不等式

    在概率论中,马尔可夫不等式给出了随机变量的非负函数大于或等于某个正常数 ϵ epsilon ϵ 的概率的上限 下图来自:Markov inequality 下图为任一分布的概率密度函数图像 图片来自:Mathematical Foundations of Computer Networking: Probability a a a 越大,阴影部分的面积越小,即概率越小 使

    2024年02月04日
    浏览(18)
  • 概论_第5章_切比雪夫不等式

    一 定理: 切比雪夫不等式 设随机变量 X的期望E(X) 及方差D(X)存在, 则对任意小正数 0, 有 P{|X-E(X)| } , 或者 P{|X-E(X)| } 1- 其中 念伊普西龙 可见, 要使用切比雪夫不等式, 先要算出E(X), D(X), 再根据已知条件算出 , 然后才能套用公式 。 例题 设电站供电网有 1万盏灯, 夜晚每一

    2024年02月12日
    浏览(16)
  • Open3D 最小二乘拟合平面(SVD分解法)

    本文由CSDN点云侠原创,原文链接。爬虫网站自重。 本文实现矩阵奇异值分解方法的最小二乘拟合平面。原理如下:   对于得到的

    2024年02月13日
    浏览(20)
  • 25、基于原型的切比雪夫低通滤波器匹配设计理论

    低通滤波器匹配网络其实就是在滤波的基础上增加了一个阻抗变换的作用,其设计参数包含阻抗变换比、设计带宽参数等等,因为其良好的匹配特性所以经常使用在功率放大器的设计之中。如MTT文章: Design of Highly Efficient Broadband Class-E Power Amplifier Using Synthesized Low-Pass Matchin

    2024年02月08日
    浏览(25)
  • 3D点云处理:用SVD分解法和最小二乘法拟合平面点云,求解平面方程

    本文主要介如何用SVD分解法和最小二乘法拟合平面点云,包含原理推导和代码 将空间中的离散点拟合为一个平面,就是使离散点到某个平面距离和最小的问题,可以将求解过程看作最优化的过程。 一个先验知识为拟合平面一定经过离散点的质心(离散点坐标的平均值)。平

    2024年02月03日
    浏览(77)
  • 24、基于原型的切比雪夫低通滤波器设计理论(插入损耗法)

    基于原型的滤波器设计是设计的一个基础的方法,虽然在现在有更加强大的自动化工具能够取代它,但是如果要进行理论研究仍需要对其有所了解。 写此文的初衷并非是介绍切比雪夫低通滤波器设计理论,而是发现国内有一些文章把 低通滤波器网络 和 低通滤波器匹配网络

    2023年04月22日
    浏览(12)
  • 3D点云平面拟合算法

    假设你有一组 3D 中的 n 个点,并且想要为它们拟合一个平面。 在本文中,我将推导出一个简单的、数值稳定的方法,并提供它的源代码。 听起来很好玩? 我们开始吧! NSDT工具推荐 : Three.js AI纹理开发包 - YOLO合成数据生成器 - GLTF/GLB在线编辑 - 3D模型格式在线转换 

    2024年02月02日
    浏览(20)
  • 深度神经网络的数学原理:基于超平面、半空间与线性区域的表示

    以前的文章主要描述了神经网络,即多层感知机、全连接模型的运行原理,还是以实验为主,数学描述为辅的方式,这篇文章以纯数学的视角来描述神经网络的运行原理,主要以前馈过程为主(反向传播的动力学过程还是比较复杂,正向过程还未完全研究清楚,暂时还未考虑

    2024年02月06日
    浏览(21)
  • OpenCV学堂 | CV开发者必须懂的9种距离度量方法,内含欧氏距离、切比雪夫距离等(建议收藏)

    本文来源公众号“ OpenCV学堂 ”,仅用于学术分享,侵权删,干货满满。 原文链接:CV开发者必须懂的9种距离度量方法,内含欧氏距离、切比雪夫距离等 在数据挖掘中,我们经常需要计算样本之间的相似度,通常的做法是计算样本之间的距离。在本文中,数据科学家 Maarten

    2024年02月20日
    浏览(24)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包