【数理知识】求刚体旋转矩阵和平移矩阵,已知 N>=3 个点在前后时刻的坐标,且这 N>=3 点间距离始终不变代表一个刚体

序号内容
1【数理知识】自由度 degree of freedom 及自由度的计算方法
2【数理知识】刚体 rigid body 及刚体的运动
3【数理知识】刚体基本运动,平动,转动
4【数理知识】向量数乘,内积,外积,matlab代码实现
5【数理知识】协方差,随机变量的的协方差,随机变量分别是单个数字和向量时的协方差
6【数理知识】已知 N>=3 个点在前后时刻的坐标,求刚体平移矩阵,旋转矩阵,且这 N>=3 点间距离始终不变代表一个刚体

文章目录

  • 1 解决流程
    • 1. 找寻质心(Centroid)
    • 2. 奇异值分解(SVD)
    • 3. 通过协方差矩阵得到旋转矩阵
    • 4. 计算平移矩阵
  • 2 举例验证 1
    • 1. 找寻质心(Centroid)
    • 2. 计算协方差矩阵
    • 3. 奇异值分解
    • 4. 计算平移矩阵
  • Ref

存在有 N ≥ 3 N\ge 3 N3 个点,它们两两之间距离始终不变,这就满足了可代表一个刚体的条件。同时,已知这 N ≥ 3 N\ge 3 N3 个点在前后时刻的坐标,如何求对应刚体的平移矩阵,旋转矩阵?

如下图所示,对应点的颜色相同, R R R 是旋转, t t t 是平移。我们希望找到能将数据集 A A A 中的点对齐到数据集 B B B 的最佳旋转和平移。这种变换有时被称为欧几里得变换(Euclidean)或刚性变换(Rigid transform),因为它保留了形状和大小。这与仿射变换不同,后者包括缩放和剪切。

这个问题尤其出现在三维点云数据注册等任务中,因为这些数据是从三维激光扫描仪或流行的 Kinect 设备等硬件中获取的。

在这里插入图片描述

接下来的描述中,我们为了和图中情况保持一致,我们都假设 N = 3 N = 3 N=3

1 解决流程

旋转和平移的方程式可以表示为如下形式:

R A + t = B \begin{aligned} RA + t = B \end{aligned} RA+t=B

最终目的是求取最合适的 R R R t t t

至于为什么可以这么表示,请参考文章开头所提及到的其他文章。

1. 找寻质心(Centroid)

这一步也比较简单,质心就是 N = 3 N = 3 N=3 个数据点的平均值

Centroid A = 1 3 ∑ k = 1 3 A k Centroid B = 1 3 ∑ k = 1 3 B k \begin{aligned} \text{Centroid}_{A} &= \frac{1}{3} \sum_{k=1}^{3} A_k \\ \text{Centroid}_{B} &= \frac{1}{3} \sum_{k=1}^{3} B_k \end{aligned} CentroidACentroidB=31k=13Ak=31k=13Bk

其中 A k A_k Ak B k B_k Bk 分别表示在数据集 A A A B B B 中第 k k k 个数据点的坐标。


2. 奇异值分解(SVD)

有几种方法可以找到点之间的最佳旋转。最简单的方法是使用奇异值分解(SVD),因为许多编程语言(Matlab、Octave、使用 LAPACK 的 C 语言、使用 OpenCV 的 C++ 语言…)都可以使用这个函数。SVD 就像线性代数中的一根神奇魔杖,可以解决各种数值问题。这里不会详细介绍它的工作原理,而会介绍如何使用它。你只需要知道,SVD 可以将一个矩阵(称作 E E E)分解/因式分解为另外 3 个矩阵,即

[ U , S , V ] = SVD ( E ) E = U S V T \begin{aligned} [U, S, V] &= \text{SVD} (E) \\ E &= U S V^\text{T} \end{aligned} [U,S,V]E=SVD(E)=USVT

如果 E E E 是方阵,那么 U 、 S U、S US V V V 的大小也相同。

3. 通过协方差矩阵得到旋转矩阵

要找到最佳旋转方式,我们首先要重新调整两个数据集的中心,使两个中心点都位于原点,如下图所示。

在这里插入图片描述

这样就去除了平移部分,只剩下旋转部分需要处理。下一步是累加一个矩阵(称为 H H H),然后使用 SVD 求出旋转,如下所示:

H = ( A − Centroid A ) ( B − Centroid B ) T [ U , S , V ] = SVD ( H ) R = V U T \begin{aligned} H &= (A - \text{Centroid}_{A})(B - \text{Centroid}_{B})^\text{T} \\ [U, S, V] &= \text{SVD} (H) \\ R &= V U^\text{T} \end{aligned} H[U,S,V]R=(ACentroidA)(BCentroidB)T=SVD(H)=VUT

其中, H H H 是我们熟悉的协方差矩阵。 A − Centroid A A - \text{Centroid}_{A} ACentroidA 是用 A A A 减去 Centroid A \text{Centroid}_{A} CentroidA 中的每一列的操作。

需要注意的一点是,要正确计算 H H H。它最终应该是一个 3 × 3 3 \times 3 3×3 矩阵,而不是一个 N × N N \times N N×N 矩阵(这里 N N N 是指点的数量,而 3 3 3 是指数据的坐标 [ x , y , z ] [x,y,z] [x,y,z] 维度是 3 3 3)。注意转置符号。它是在两个矩阵之间进行乘法运算,这两个矩阵的实际维数分别是 3 × N 3 \times N 3×N N × 3 N \times 3 N×3。乘法的顺序也很重要,如果换一种方法,就会变成是从 B B B A A A 的旋转。

4. 计算平移矩阵

得到旋转矩阵 R R R 后,平移矩阵 t t t 也就变得简单了。把质心代入开篇咱们提到的方程,那么有

R A + t = B R × Centroid A + t = Centroid B t = Centroid B − R × Centroid A \begin{aligned} RA + t &= B \\ R \times \text{Centroid}_A + t &= \text{Centroid}_B \\ t &= \text{Centroid}_B - R \times \text{Centroid}_A \end{aligned} RA+tR×CentroidA+tt=B=CentroidB=CentroidBR×CentroidA


2 举例验证 1

假设有 3 3 3 个相对位置保持不变的点,已知它们在数据集合 A A A 和数据集合 B B B 中的位置,然后计算旋转矩阵 R R R 和平动矩阵 t t t

在数据集合 A A A 中:
1 1 1 的位置为: A 1 = ( 1 , 2 , 3 ) A_1 = (1, 2, 3) A1=(1,2,3)
2 2 2 的位置为: A 2 = ( 4 , 5 , 6 ) A_2 = (4, 5, 6) A2=(4,5,6)
3 3 3 的位置为: A 3 = ( 7 , 8 , 9 ) A_3 = (7, 8, 9) A3=(7,8,9)

在数据集合 B B B 中:
1 1 1 的位置为: B 1 = ( 2 , 3 , 4 ) B_1 = (2, 3, 4) B1=(2,3,4)
2 2 2 的位置为: B 2 = ( 5 , 6 , 7 ) B_2 = (5, 6, 7) B2=(5,6,7)
3 3 3 的位置为: B 3 = ( 8 , 9 , 10 ) B_3 = (8, 9, 10) B3=(8,9,10)

计算思路为:

  • 先计算平移:通过求取这些点在两个时刻的质心位置,然后求差来得到平移矩阵

1. 找寻质心(Centroid)

这一步也比较简单,直接代入样本数据

Centroid A = ( 1 + 4 + 7 , 2 + 5 + 8 , 3 + 6 + 9 ) 3 = ( 1 + 4 + 7 3 , 2 + 5 + 8 3 , 3 + 6 + 9 3 ) = ( 4 , 5 , 6 ) Centroid B = ( 2 + 5 + 8 , 3 + 6 + 9 , 4 + 7 + 10 ) 3 = ( 2 + 5 + 8 3 , 3 + 6 + 9 3 , 4 + 7 + 10 3 ) = ( 5 , 6 , 7 ) \begin{aligned} \text{Centroid}_{A} &= \frac{(1+4+7, 2+5+8, 3+6+9)}{3} = (\frac{1 + 4 + 7}{3}, \frac{2 + 5 + 8}{3}, \frac{3 + 6 + 9}{3}) = (4, 5, 6) \\ \text{Centroid}_{B} &= \frac{(2+5+8, 3+6+9, 4+7+10)}{3} = (\frac{2 + 5 + 8}{3}, \frac{3 + 6 + 9}{3}, \frac{4 + 7 + 10}{3}) = (5, 6, 7) \end{aligned} CentroidACentroidB=3(1+4+7,2+5+8,3+6+9)=(31+4+7,32+5+8,33+6+9)=(4,5,6)=3(2+5+8,3+6+9,4+7+10)=(32+5+8,33+6+9,34+7+10)=(5,6,7)

2. 计算协方差矩阵

根据公式

Cov ( X , Y ) i j = ∑ k n = 3 ( x k i − x ˉ i ) ( y k j − y ˉ i ) n − 1 \begin{aligned} \text{Cov} (X,Y)_{ij} &= \frac{\sum_k^{n=3} (x_{ki} - \bar{x}_i)(y_{kj} - \bar{y}_i)}{n-1} \end{aligned} Cov(X,Y)ij=n1kn=3(xkixˉi)(ykjyˉi)

可以得到协方差矩阵

Cov ( X , Y ) = [ 3 3 3 3 3 3 3 3 3 ] \begin{aligned} \text{Cov} (X,Y) &= \left[\begin{matrix} 3 & 3 & 3 \\ 3 & 3 & 3 \\ 3 & 3 & 3 \\ \end{matrix}\right] \end{aligned} Cov(X,Y)= 333333333

关于协方差矩阵的原理和求解方法,可参考文章:【数理知识】协方差,随机变量的的协方差,随机变量分别是单个数字和向量时的协方差。

3. 奇异值分解

可以直接使用 Matlab 进行奇异值分解,可以得到

U = [ − 0.5774 0.8165 − 0.0000 − 0.5774 − 0.4082 − 0.7071 − 0.5774 − 0.4082 0.7071 ] , S = [ 9 0 0 0 0 0 0 0 0 ] , V = [ − 0.5774 0.8165 0 − 0.5774 − 0.4082 − 0.7071 − 0.5774 − 0.4082 0.7071 ] U = \left[\begin{matrix} -0.5774 & 0.8165 & -0.0000 \\ -0.5774 & -0.4082 & -0.7071 \\ -0.5774 & -0.4082 & 0.7071 \\ \end{matrix}\right], S = \left[\begin{matrix} 9 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{matrix}\right], V = \left[\begin{matrix} -0.5774 & 0.8165 & 0 \\ -0.5774 & -0.4082 & -0.7071 \\ -0.5774 & -0.4082 & 0.7071 \\ \end{matrix}\right] U= 0.57740.57740.57740.81650.40820.40820.00000.70710.7071 ,S= 900000000 ,V= 0.57740.57740.57740.81650.40820.408200.70710.7071

R = V U T = [ 1 0 0 0 1 0 0 0 1 ] R = V U^\text{T} = \left[\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{matrix}\right] R=VUT= 100010001

至此得到了旋转矩阵。

4. 计算平移矩阵

t = Centroid B − R × Centroid A = ( 1 , 1 , 1 ) \begin{aligned} t &= \text{Centroid}_B - R \times \text{Centroid}_A &= (1, 1, 1) \end{aligned} t=CentroidBR×CentroidA=(1,1,1)


Ref

  1. FINDING OPTIMAL ROTATION AND TRANSLATION BETWEEN CORRESPONDING 3D POINTS
  2. 从3组对应点中求得最佳的旋转和平移变换

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/61257.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

揭秘程序员最喜欢的5个高薪工作

大家好,这里是程序员晚枫。想了解更多精彩内容,快来关注程序员晚枫 今天给大家推荐5个适合程序员的高薪岗位。 01 推荐岗位 以下是5个工资最高的程序员工作: 数据科学家:数据科学家是负责数据收集、处理、分析和报告的专业人员。…

点击编辑变完成

<template><div><button click"dialogshowtrue">添加部门</button><div>部门列表</div><el-table ref"multipleTable" :data"form" tooltip-effect"dark" style"width: 100%">&l…

程序环境和预处理(含C语言程序的编译+链接)--1

&#x1f389;个人名片&#xff1a; &#x1f43c;作者简介&#xff1a;一名乐于分享在学习道路上收获的大二在校生 &#x1f43b;‍❄个人主页&#x1f389;&#xff1a;GOTXX &#x1f43c;个人WeChat&#xff1a;ILXOXVJE &#x1f43c;本文由GOTXX原创&#xff0c;首发CSDN…

【Python】数据分析+数据挖掘——掌握Python和Pandas中的单元格替换操作

1. 前言 数据处理和清洗是数据分析和机器学习中至关重要的步骤。在数据处理过程中&#xff0c;我们经常需要对数据集进行清洗和转换&#xff0c;其中单元格替换是一个常用的技术。Python作为一种功能强大且灵活的编程语言&#xff0c;为数据处理提供了丰富的工具和库。Pandas库…

tinkerCAD案例:31. 3D 基元形状简介

tinkerCAD案例&#xff1a;31. 3D 基元形状简介 1 将一个想法从头脑带到现实世界是一次令人兴奋的冒险。在 Tinkercad 中&#xff0c;这将从一个新的设计开始。 在新设计中&#xff0c;简单的原始形状可以通过不同的方式组合成更复杂的形状。 在这个项目中&#xff0c;你将探索…

网络安全设备及部署

什么是等保定级&#xff1f; 之前了解了下等保定级&#xff0c;接下里做更加深入的探讨 文章目录 一、网路安全大事件1.1 震网病毒1.2 海康威视弱口令1.3 物联网Mirai病毒1.4 专网 黑天安 事件1.5 乌克兰停电1.6 委内瑞拉电网1.7 棱镜门事件1.8 熊猫烧香 二、法律法规解读三、安…

C语言基础知识点一

C语言基础知识点一&#xff1a; 1.数据类型 2.bool类型&#xff1a; 使用bool时时&#xff0c;需要增加<stdbool.h>头文件。 说明&#xff1a;bool 类型只有非零&#xff08;true&#xff09;和零&#xff08;false&#xff09;两种值。 如: if&#xff08;-1&#xf…

【供电并接电路】2021-12-31

缘由在实际应用中mos管的导通问题 - 电路设计论坛 - 电子技术论坛 - 广受欢迎的专业电子论坛!MOS管实际应用中的开通问题-硬件开发-CSDN问答 设计的思想是&#xff0c;当CH_5.2V有5.2V电压时&#xff0c;Q7导通&#xff0c;Q6产生UGS​压差&#xff0c;从而使Q6导通&#xff0c…

SpringBoot+AOP+Redission实战分布式锁

文章目录 前言一、Redission是什么&#xff1f;二、使用场景三、代码实战1.项目结构2.类图3.maven依赖4.yml5.config6.annotation7.aop8.model9.service 四、单元测试总结 前言 在集群环境下非单体应用存在的问题&#xff1a;JVM锁只能控制本地资源的访问&#xff0c;无法控制…

Python中的PDF文本提取:使用fitz和wxPython库(带进度条)

引言&#xff1a; 处理大量PDF文档的文本提取任务可能是一项繁琐的工作。本文将介绍一个使用Python编写的工具&#xff0c;可通过简单的操作一键提取大量PDF文档中的文本内容&#xff0c;极大地提高工作效率。 import wx import pathlib import fitzclass PDFExtractor(wx.Fr…

2.运行Python

在完成上一篇的环境安装后,本篇文章我们运行python及简单的命令使用 在安装之后我们可以在最近安装的程序中点击快捷方式启动python 第一种方式 python自带的开发环境 双击安装好后的应用程序 第二种方式 命令行 将print(“hello word”) 写入到文件,文件后缀名以.py结…

Windows环境下VSCode安装PlatformIO Cero报错ERROR: HTTP error 403 while getting

安装PlatformIO插件成功&#xff0c;初始化失败 错误信息判断问题尝试访问https://pypi.tuna.tsinghua.edu.cn/simple/platformio/成功点击文件后报错如下&#xff1a; 解决问题- 换源 &#xff08; Windows下有两个地方需要更改&#xff09;cmd命令行Pip文件 总结&#xff1a;…

Android Gradle 骚操作,将两个项目合并到一个项目中

1. 前言 在工作中&#xff0c;由于各种原因&#xff0c;导致需要将两个可单独运行的App项目&#xff0c;合并到一个git仓库里&#xff0c;且单独的App项目里还有其他Module模块。 如果只是将两个项目复制到同一个文件夹下&#xff0c;还是得单独打开各个项目&#xff0c;是很不…

python自动化运维常用模块,python自动化运维项目

大家好&#xff0c;本文将围绕python自动化运维需要掌握的技能展开说明&#xff0c;python自动化运维快速入门 pdf是一个很多人都想弄明白的事情&#xff0c;想搞清楚python自动化运维常用模块需要先了解以下几个事情。 这篇文章主要介绍了一个有趣的事情&#xff0c;具有一定借…

b 树和 b+树的理解

项目场景&#xff1a; 图灵奖获得者&#xff08;Niklaus Wirth &#xff09;说过&#xff1a; 程序 数据结构 算法&#xff0c; 也就说我们无时无刻 都在和数据结构打交道。 只是作为 Java 开发&#xff0c;由于技术体系的成熟度较高&#xff0c;使得大部分人认为&#xff1…

英特尔傲腾CAS报错unknown error cache acceleration software could not start cache

英特尔傲腾CAS报错unknown error cache acceleration software could not start cache 文章目录 英特尔傲腾CAS报错unknown error cache acceleration software could not start cache我是怎么遇到这个问题的我是如何解决的实验步骤打Primo Cache蓝屏补丁拔掉原来的系统盘开关机…

WSL2安装CentOS7和CentOS8

提示&#xff1a;文章写完后&#xff0c;目录可以自动生成&#xff0c;如何生成可参考右边的帮助文档 文章目录 前言一、下载ZIP包&#xff1f;二、安装1.打开Windows子系统支持2.安装到指定位置3.管理虚拟机4.配置虚拟机1.配置国内源2.安装软件3.安装第三方源 5.配置用户1.创建…

嵌入式开发学习(STC51-12-I2C/IIC)

内容 在数码管右3位显示数字&#xff0c;从0开始&#xff0c;按K1键将数据写入到EEPROM内保存&#xff0c;按K2键读取EEPROM内保存的数据&#xff0c;按K3键显示数据加1&#xff0c;按K4键显示数据清零&#xff0c;最大能写入的数据是255&#xff1b; I2C介绍 I2C简介 I2C&…

windows下的txt文档,传到ubuntu后,每行后面出现^M,怎么处理?

问题背景&#xff1a;windows下pycharm生成的txt文档&#xff0c;传到ubuntu后&#xff0c;每行后面出现^M 用vim打开显示 使用cat -A filename显示如下 参考https://www.lmlphp.com/user/16697/article/item/579325/给出的几种方法 方法一、dos2unix filename。服务器没装…

esp32c3 xiao 脚本记录

oled显示网络时间, wifi链接网络 // ntp_get_date.h #include "time.h"String week[8] {"Sun", "Mon", "Tues", "Wednes", "Thur", "Fri", "Sat" };void printLocalTime(Adafruit_SSD1306 …