《点云处理》平面拟合

前言

在众多点云处理算法中,其中关于平面拟合的算法十分广泛。本篇内容主要是希望总结归纳各类点云平面拟合算法,并且将代码进行梳理保存。

环境:

VS2019 + PCL1.11.1

1.RANSAC

使用ransac对平面进行拟合是非常常见的用法,PCL库中就有RANSAC拟合平面的实现代码,而且还集成了 两种拟合平面的代码。
方法一:

/// <summary>
/// 使用PCL中集成的RANSAC方法进行平面拟合
/// </summary>
/// <param name="cloud_in">输入待拟合的点云</param>
/// <param name="inliers">RANSAC拟合得到的内点</param>
/// <param name="coefficients">得到的平面方程参数</param>
/// <param name="iterations">平面拟合最大迭代次数</param>
/// <param name="threshold">RANSAC拟合算法距离阈值</param>
void RANSAC(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud_in, std::vector<int>& inliers, Eigen::VectorXf& coefficients,const int& iterations, const double& threshold)
{inliers.clear();                                          // 用于存放内点索引的vectorpcl::shared_ptr<pcl::SampleConsensusModelPlane<pcl::PointXYZ>> model(new pcl::SampleConsensusModelPlane<pcl::PointXYZ>(cloud_in)); //定义待拟合平面的model,并使用待拟合点云初始化pcl::RandomSampleConsensus<pcl::PointXYZ> ransac(model);  // 定义RANSAC算法模型ransac.setDistanceThreshold(threshold);                   // 设定阈值ransac.setMaxIterations(iterations);                      // 设置最大迭代次数ransac.setNumberOfThreads(10);                            // 设置线程数量ransac.computeModel();                                    // 拟合 ransac.getInliers(inliers);                               // 获取内点索引Eigen::VectorXf params;                                   // 第一次得到的平面方程ransac.getModelCoefficients(params);                      // 获取拟合平面参数,对于平面ax+by_cz_d=0,params分别按顺序保存a,b,c,dmodel->optimizeModelCoefficients(inliers, params, coefficients);      // 优化平面方程参数std::vector<double> vDistance;                            // 用于存储每个点到拟合平面的距离的vector容器model->getDistancesToModel(coefficients, vDistance);      // 得到每个点到平面的距离的集合std::cout << "params: " << params[0] << ", " << params[1] << ", " << params[2] << ", " << params[3] << std::endl;std::cout << "coefficients: " << coefficients[0] << ", " << coefficients[1] << ", " << coefficients[2] << ", " << coefficients[3] << std::endl;
}

上述代码需要注意的是,一定需要model->optimizeModelCoefficients(inliers, params, coefficients); 通过这一句代码去优化平面参数。优化前后差别很大,如下图所示:
在这里插入图片描述
方法二:

/// <summary>
/// 使用PCL中集成的用于点云分割的RANSAC方法进行平面拟合
/// </summary>
/// <param name="cloud_in">输入待拟合的点云</param>
/// <param name="inliers">RANSAC拟合得到的内点</param>
/// <param name="coefficients">得到的平面方程参数</param>
/// <param name="iterations">平面拟合最大迭代次数</param>
/// <param name="threshold">RANSAC拟合算法距离阈值</param>
void SEG_RANSAC(const pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud_in, pcl::PointIndices::Ptr& inliers, Eigen::VectorXf& coefficients,const int& iterations, const double& threshold)
{if (inliers == nullptr) inliers.reset(new pcl::PointIndices);pcl::ModelCoefficients::Ptr coefficients_m(new pcl::ModelCoefficients);pcl::SACSegmentation<pcl::PointXYZ> seg;seg.setOptimizeCoefficients(true);seg.setModelType(pcl::SACMODEL_PLANE);seg.setMethodType(pcl::SAC_RANSAC);seg.setMaxIterations(iterations);                     // 设置最大迭代次数seg.setDistanceThreshold(threshold);                  // 设定阈值seg.setNumberOfThreads(10);                           // 设置线程数量seg.setInputCloud(cloud_in);seg.segment(*inliers, *coefficients_m);coefficients.resize(4);coefficients[0] = coefficients_m->values[0]; coefficients[1] = coefficients_m->values[1];coefficients[2] = coefficients_m->values[2]; coefficients[3] = coefficients_m->values[3];std::cout << "SEG coefficients: " << coefficients[0] << ", " << coefficients[1] << ", " << coefficients[2] << ", " << coefficients[3] << std::endl;
}

方法二其实也提到了优化平面参数,seg.setOptimizeCoefficients(true);这句代码就是比较关键的,需要加上。
上述两种方法的运行结果如下,从结果和运行时间来看,两种方法几乎一致。
在这里插入图片描述
两种方法中都有一个设置线程的功能setNumberOfThreads(10); 但是这句话根本没有起到作用,下图中有提示[pcl::RandomSampleConsensus::computeModel] Parallelization is requested, but OpenMP 3.1 is not available! Continuing without parallelization.其实在VS中已经开启了openmp,也包含了头文件#include <omp.h>,只不过版本不对应,该问题也查阅了很久,没查到解决办法。

VS中开启openmp加速的配置方法:
在这里插入图片描述

2.最小二乘法拟合平面

除了RANSAC之外,另一种十分常见的拟合方程的方法是最小二乘法,这里就不过多介绍相关的理论知识,这里提到了相关的代码和原理

当然还是要贴上封装的C++代码实现:

/// <summary>
/// 使用最小二乘法拟合平面:Ax + By + Cz + D = 0 */   //拉格朗日乘子法 https://zhuanlan.zhihu.com/p/390294059
/// </summary>
/// <param name="xjData">输入待拟合平面的点云</param>
/// <param name="coefficients">输出拟合的平面方程</param>
/// <returns>输入点云点数符合要求,返回true,否则返回false</returns>
bool LeastSquare(const pcl::shared_ptr<pcl::PointCloud<pcl::PointXYZ>>& xjData, Eigen::VectorXf& coefficients)
{int count = xjData->points.size();if (count < 3) return false;double meanX = 0, meanY = 0, meanZ = 0;double meanXX = 0, meanYY = 0, meanZZ = 0;double meanXY = 0, meanXZ = 0, meanYZ = 0;for (int i = 0; i < count; i++){meanX += xjData->points[i].x;meanY += xjData->points[i].y;meanZ += xjData->points[i].z,meanXX += xjData->points[i].x * xjData->points[i].x;meanYY += xjData->points[i].y * xjData->points[i].y;meanZZ += xjData->points[i].z * xjData->points[i].z;meanXY += xjData->points[i].x * xjData->points[i].y;meanXZ += xjData->points[i].x * xjData->points[i].z;meanYZ += xjData->points[i].y * xjData->points[i].z;}meanX /= count; meanY /= count; meanZ /= count;meanXX /= count; meanYY /= count; meanZZ /= count;meanXY /= count; meanXZ /= count; meanYZ /= count;/* eigenvector */Eigen::Matrix3d eMat;eMat(0, 0) = meanXX - meanX * meanX; eMat(0, 1) = meanXY - meanX * meanY; eMat(0, 2) = meanXZ - meanX * meanZ;eMat(1, 0) = meanXY - meanX * meanY; eMat(1, 1) = meanYY - meanY * meanY; eMat(1, 2) = meanYZ - meanY * meanZ;eMat(2, 0) = meanXZ - meanX * meanZ; eMat(2, 1) = meanYZ - meanY * meanZ; eMat(2, 2) = meanZZ - meanZ * meanZ;Eigen::EigenSolver<Eigen::Matrix3d> xjMat(eMat);           // 求取矩阵特征值和特征向量的函数EigenSolverEigen::Matrix3d eValue = xjMat.pseudoEigenvalueMatrix();   // 获取矩阵伪特征值 3*3Eigen::Matrix3d eVector = xjMat.pseudoEigenvectors();      // 获取矩阵伪特征向量 3*1/* the eigenvector corresponding to the minimum eigenvalue */double v1 = eValue(0, 0); double v2 = eValue(1, 1); double v3 = eValue(2, 2);int minNumber = 0;if ((abs(v2) <= abs(v1)) && (abs(v2) <= abs(v3))){minNumber = 1;}if ((abs(v3) <= abs(v1)) && (abs(v3) <= abs(v2))){minNumber = 2;}double A = eVector(0, minNumber); double B = eVector(1, minNumber); double C = eVector(2, minNumber);double length = sqrt(A * A + B * B + C * C);A /= length; B /= length; C /= length;double D = -(A * meanX + B * meanY + C * meanZ);/* result */if (C < 0){A *= -1.0; B *= -1.0; C *= -1.0; D *= -1.0;}coefficients.resize(4);coefficients[0] = A; coefficients[1] = B; coefficients[2] = C; coefficients[3] = D;std::cout << "LS coefficients: " << coefficients[0] << ", " << coefficients[1] << ", " << coefficients[2] << ", " << coefficients[3] << std::endl;
}

同一份点云,使用最小二乘法拟合结果如下:
在这里插入图片描述
最小二乘法拟合平面方程和上述RANSAC拟合结果还是有些差别的,但是RANSAC运行时间约为0.02s,而最小二乘法运行时间需要0.0013s,最小二乘法运行时间要比RANSAC快很多。如果针对一份没有太多噪点,很平滑干净的点云,可以使用最小二乘法去拟合,但是如果是有噪点的点云,追求准确率的情况下,则需要使用RANSAC进行拟合。

3.SVD分解的方法求平面方程

可以通过Eigen实现使用SVD分解的方法进行平面方程求解

/// <summary>
/// 通过SVD分解的方法求拟合平面
/// </summary>
/// <param name="pcl_cloud">输入待拟合平面点云</param>
/// <param name="coefficients">输出拟合平面的参数</param>
void PlaneEigenSVD(const pcl::PointCloud<pcl::PointXYZ>& pcl_cloud, Eigen::VectorXf& coefficients)
{if (pcl_cloud.points.size() < 3) return;// Convert PCL point cloud to Eigen matrixEigen::Matrix<float, 3, Eigen::Dynamic> coord(3, pcl_cloud.size());for (size_t i = 0; i < pcl_cloud.size(); ++i){coord.col(i) << pcl_cloud[i].x, pcl_cloud[i].y, pcl_cloud[i].z;}// Calculate centroidEigen::Vector4f centroid;pcl::compute3DCentroid(pcl_cloud, centroid);// Subtract centroidcoord.row(0).array() -= centroid(0);coord.row(1).array() -= centroid(1);coord.row(2).array() -= centroid(2);// Perform singular value decomposition (SVD)Eigen::JacobiSVD<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic>> svd(coord, Eigen::ComputeThinU | Eigen::ComputeThinV);Eigen::Vector3f plane_normal = svd.matrixU().rightCols<1>();// Create PCL ModelCoefficientscoefficients.resize(4);coefficients[0] = plane_normal(0);coefficients[1] = plane_normal(1);coefficients[2] = plane_normal(2);coefficients[3] = -(plane_normal(0) * centroid(0) + plane_normal(1) * centroid(1) + plane_normal(2) * centroid(2));if (coefficients[2] < 0) coefficients = -coefficients;std::cout << "SVD coefficients: " << coefficients[0] << ", " << coefficients[1] << ", " << coefficients[2] << ", " << coefficients[3] << std::endl;return;
}

运行结果如下:
在这里插入图片描述
得到的结果其实和最小二乘法拟合得到的结果一样,但是耗时更长一些。

4.霍夫变换进行平面拟合

其实使用霍夫变换去拟合几何模型也是一件非常常见的方法,但是目前没有相关代码进行测试验证,如果后续有找到实现方法也会更新在这里。

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

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

医疗智能化革命:AI技术引领医疗领域的创新进程

一、“AI”医疗的崛起 随着人工智能&#xff08;AI&#xff09;技术的崛起&#xff0c;"AI"医疗正在以惊人的速度改变着医疗行业的面貌。AI作为一种强大的工具&#xff0c;正在为医疗领域带来前所未有的创新和突破。它不仅在医学影像诊断、病理学分析和基因组学研究等…

设计模式——状态模式

引言 状态模式是一种行为设计模式&#xff0c; 让你能在一个对象的内部状态变化时改变其行为&#xff0c; 使其看上去就像改变了自身所属的类一样。 问题 状态模式与有限状态机 的概念紧密相关。 其主要思想是程序在任意时刻仅可处于几种有限的状态中。 在任何一个特定状态中…

手拉手EasyExcel极简实现web上传下载(全栈)

环境介绍 技术栈 springbootmybatis-plusmysqleasyexcel 软件 版本 mysql 8 IDEA IntelliJ IDEA 2022.2.1 JDK 1.8 Spring Boot 2.7.13 mybatis-plus 3.5.3.2 EasyExcel是一个基于Java的、快速、简洁、解决大文件内存溢出的Excel处理工具。 他能让你在不用考虑性…

CentOS7安装教程

1.准备工作 安装虚拟机软件 VMware Workstation 17 Player下载CentOS7 镜像 2. 安装VMware Workstation 17 Player 官网 下载链接 下载好了安装包&#xff0c;双击安装包&#xff0c;傻瓜式安装一直下一步&#xff0c;安装即可。 3. 安装CentOS7 官网 推荐下载地址&…

Tekton 构建容器镜像

Tekton 构建容器镜像 介绍如何使用 Tektonhub 官方 kaniko task 构建docker镜像&#xff0c;并推送到远程dockerhub镜像仓库。 kaniko task yaml文件下载地址&#xff1a;https://hub.tekton.dev/tekton/task/kaniko 查看kaniko task yaml内容&#xff1a; 点击Install&…

RabbitMQ 消息持久化

默认情况下&#xff0c;exchange、queue、message 等数据都是存储在内存中的&#xff0c;这意味着如果 RabbitMQ 重启、关闭、宕机时所有的信息都将丢失。 RabbitMQ 提供了持久化来解决这个问题&#xff0c;持久化后&#xff0c;如果 RabbitMQ 发送 重启、关闭、宕机&#xff…

05 动态渲染数据

概述 实际上动态渲染数据&#xff0c;在《使用CDN开发Vue3项目》中就已经学习过了&#xff0c;核心代码如下&#xff1a; <div id"vue-app">{{text}}</div> <script src"https://cdn.staticfile.org/vue/3.0.5/vue.global.js"></sc…

【运维笔记】Hyperf正常情况下Xdebug报错死循环解决办法

问题描述 在使用hyperf进行数据库迁移时&#xff0c;迁移报错&#xff1a; 查看报错信息&#xff0c;错误描述是Xdebug检测到死循环&#xff0c;可是打印的堆栈确实正常堆栈&#xff0c;没看到死循环。 寻求解决 gpt 说的跟没说一样。。 google一下 直接把报错信息粘贴上去…

I.MX RT1170双核学习(3):多核管理之MCMGR源码分析详解

本文通过SDK中最简单的hello_world例程来说明一下双核程序如何运行。在CM7和CM4的工程中都有一个MCMGR(Multicore Manager)文件夹&#xff0c;它是用来管理多核之间的操作的&#xff0c;当然也包括我们前面提到的那些寄存器的设置。 文章目录 1 MCMGR_EarlyInit1.1 MCMGR_Trigg…

MFC 程序执行流程

目录 MFC 程序启动 MFC 入口函数 程序执行流程总结 在Win32课程中WinMain由程序员自己实现&#xff0c;那么流程是程序员安排&#xff0c;但到了MFC中&#xff0c;由于MFC库实现WinMain&#xff0c;也就意味着MFC负责安排程序的流程。 MFC 程序启动 程序的启动&#xff0c;…

PCB设计规则中的经验公式_笔记

PCB设计规则中的经验公式 规则1 - 临界长度规则2 - 信号带宽与上升时间规则3- 时钟信号带宽规则4-信号传输速度规则5- 集肤 (效应) 深度规则6 - 50Ω传输线电容规则7 - 50Ω传输线电感规则8 - 回流路径电感规则9 - 地弹噪声规则10- 串行传输比特率与信号带宽规则11- PCB走线直流…

HIVE窗口函数

什么是窗口函数 hive中开窗函数通过over关键字声明&#xff1b;窗口函数&#xff0c;准确地说&#xff0c;函数在窗口中的应用&#xff1b;比如sum函数不仅可在group by后聚合&#xff0c;在可在窗口中应用&#xff1b; hive中groupby算子和开窗over&#xff0c;shuffle的逻辑…

面试 Java 算法高频题五问五答第一期

面试 Java 算法高频题五问五答第一期 作者&#xff1a;程序员小白条&#xff0c;个人博客 相信看了本文后&#xff0c;对你的面试是有一定帮助的&#xff01; ⭐点赞⭐收藏⭐不迷路&#xff01;⭐ 1&#xff09;括号生成: 数字 n 代表生成括号的对数&#xff0c;请你设计一个…

OpenSergo Dubbo 微服务治理最佳实践

*作者&#xff1a;何家欢&#xff0c;阿里云 MSE 研发工程师 Why 微服务治理&#xff1f; 现代的微服务架构里&#xff0c;我们通过将系统分解成一系列的服务并通过远程过程调用联接在一起&#xff0c;在带来一些优势的同时也为我们带来了一些挑战。 如上图所示&#xff0c;可…

C语言之枚举类型

目录 枚举类型 枚举常量 枚举类型的特征 命名空间 本节我们来学习表示一定整数值的集合的枚举类型。 枚举类型 老样子&#xff0c;我们先用一段程序引出&#xff1a; /*显示所选动物的叫声*/ #include<stdio.h>enum animal {Dog, Cat, Monkey, Invalid}; /*显示狗叫…

Zotero攻略

给大家分享一下我对于Zotero的使用。 1、下载链接 Zotero | Your personal research assistant 进入后直接下载即可 2、一些好用的插件 &#xff08;1&#xff09;Zotero Connector 下载地址&#xff1a;Zotero | Connectors 超级好用&#xff01;不用一篇一篇下PDF了&am…

Redis设计与实现之事务

一、事务 Redis 通过 MULTI 、DISCARD 、EXEC 和 WATCH 四个命令来实现事务功能&#xff0c;本章首先讨 论使用 MULTI 、DISCARD 和 EXEC 三个命令实现的一般事务&#xff0c;然后再来讨论带有 WATCH 的事务的实现。 因为事务的安全性也非常重要&#xff0c;所以本章最后通过…

18个非技术面试题

请你自我介绍一下你自己&#xff1f; 这道面试题是大家在以后面试过程中会常被问到的&#xff0c;那么我们被问到之后&#xff0c;该如果回答呢&#xff1f;是说姓名&#xff1f;年龄&#xff1f;还是其他什么&#xff1f; 最佳回答提示&#xff1a; 一般人回答这个问题往往会…

neuq-acm预备队训练week 9 P1119 灾后重建

解题思路 本题可以用最短路算法——Floyd AC代码 #include<bits/stdc.h> #define inf 1e9 using namespace std; const int N 2e2 50; int n, m, q, now 0, a, b, c, t[N], G[N][N];int main() {scanf("%d%d", &n, &m);for(int i 0;i<n;i)sc…

2023新时代中国模特大赛总决赛在京落幕

12月16日&#xff0c;备受瞩目的2023新时代中国模特大赛圆满落幕。本次大赛旨在挖掘和培养具有新时代特色的模特人才&#xff0c;推动中国时尚产业的创新发展。 作为中国时尚界的重要赛事&#xff0c;新时代中国模特大赛吸引了来自全国各地的优秀模特选手45名参加全国总决赛。在…