二次曲面拟合计算点云法曲率、主曲率以及平均曲率（手写版）

2024-02-12 01:30
文章标签 计算 平均 拟合 手写 曲率 二次曲面 点云法

本文主要是介绍二次曲面拟合计算点云法曲率、主曲率以及平均曲率（手写版），希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

目录

一、相关原理及介绍

二、计算流程

三、实现代码

四、最终结果

参考文献：[1]周煜,张万兵,杜发荣等.散乱点云数据的曲率精简算法[J].北京理工大学学报,2010,30(07):785-789.DOI:10.15918/j.tbit1001-0645.2010.07.021.

                  [2]蔡志敏,王晏民,黄明.基于KD树散乱点云数据的Guass平均曲率精简算法[J].测绘通报,2013(S1):44-46.

一、相关原理及介绍

        从微分几何的角度来看，曲面曲率较高的区域，应该有较多的采样点表示，相反则应该由较少的采样点表示，曲率反映了曲面的基本特性。点云曲率特征是指在三维空间中，对于点云数据集中的每个点，通过计算其周围的几何形状来确定该点的曲率特征，曲率特征是点云数据集中非常重要的特征之一，它可以应用于点云中的分类、分割、配准、特征提取等应用中。

        曲率特征可以分为平均曲率、高斯曲率以及主曲率。平均曲率作为微分几何中一个“外在的”弯曲测量标准，对一个曲面嵌入周围空间（比如二维曲面嵌入三维欧式空间）的曲率进行了区部描述。高斯曲率作为一个描述曲面的凹凸性质的量，当这个量变化程度较大的时候表面曲面内部变化也比较大，这就表明曲面的光滑程度越低。点云中任意一点都存在某曲面 f=r(x,y,z) 逼近该点的邻域点云，一点处的曲率可用该点及其邻域点拟合的局部曲面曲率来表征。通过最小二乘拟合，可以用二次曲面来表征局部区域，计算每点处的主曲率 k{_{1}}、k_{2}，平均曲率H 及高斯曲率 K。

二、计算流程

   设点云中某一点P_{i}的K邻域点集为Q，通过邻域点建立二次曲面参数方程：

r(u,v)=\sum_{i,j=0}^{2}q_{i,j}u^{i}v^{j}                (1)

       令：

   a=[a_{00} a_{01} a_{02} a_{10} a_{11} a_{12} a_{20} a_{21} a_{22}]     b=[b_{00} b_{01} b_{02} b_{10} b_{11} b_{12} b_{20} b_{21} b_{22}]   

                        c=[c_{00} c_{01} c_{02} c_{10} c_{11} c_{12} c_{20} c_{21} c_{22}]

w=[{u^{0}}{v^{0}}\, \, {u^{0}}{v^{1}} \, \, {u^{0}}{v^{2}} \, \, {u^{1}}{v^{0}} \, \, {u^{1}}{v^{1}} \, \, {u^{1}}{v^{2}} \, \, {u^{2}}{v^{0}} \, \, {u^{2}}{v^{1}}\, \, {u^{2}}{v^{2}}]，Q=[a,b,c]        (2)

式（1）可改写成：r(u,v)=[x(u,v) \, \, y(u,v)\, \, z(u,v)]=[w^{T}a \, \,w^{T}b \, \,w^{T}c]     (3)

        引入两个矩阵：

B=\begin{bmatrix} x_{0} & y_{0} &z_{0} \\ x_{1}&y_{1} &z_{1} \\ ... &... &... \\ x_{k}&y_{k} & z_{k} \end{bmatrix}，W=\begin{bmatrix} W_{0}^{T} \\ W_{1}^{T} \\ ... \\ W_{k}^{T} \end{bmatrix}                                                                         (4)

测点P_{i}的坐标为P(x_{0}\, \, y_{0}\, \,z_{0})，选择邻域点数大于8，则逼近曲面和数据点的误差函数为：

\Omega =WQ-B                                                                                                        （5）

运用最小二乘法，推到出系数矩阵Q为：

Q=(W^{T}W)^{-1}W^{T}B                                                                                            （6）

从而求出曲面r(u,v)的偏导数r_{u}\; r_{v}\; r_{uu}\; r_{uv}\; r_{vv}，曲面的单位法矢n为：

n=\frac{r_{u} \times r_{v}}{|r_{u} \times r_{v}|}                                                                                                                    (7)

曲面的第一、第二基本量为：

E= r_{u}^{2},\; F=r_{u} r_{v},\; G=r^{2}

L= nr_{uu},\; M=nr_{uv},\; N=nr_{vv}                                                                             (8)

则，

平均曲率为：

H=\frac{EN-2FM+GL}{2(EG-F^{2})}

高斯曲率为：

K=\frac{LN-M^{2}}{EG-F^{2}}

主曲率为：

k_{1},k_{2}=H\pm \sqrt{H^{2}-K}

三、实现代码

#include <iostream>
#include <pcl/io/pcd_io.h> 
#include <pcl/point_types.h>
#include <pcl/kdtree/kdtree_flann.h>
#include <pcl/features/boundary.h>
#include <pcl/filters/normal_space.h>
#include <pcl/features/normal_3d.h>
#include <pcl/filters/radius_outlier_removal.h>
#include <vector>
#include <math.h>
#include <pcl/common/eigen.h>
#include <pcl/point_cloud.h>
#include <pcl/features/principal_curvatures.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <boost/thread/thread.hpp>using namespace std;
using namespace Eigen;//曲面拟合计算点云曲率
void surfacefitting(pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, pcl::PointCloud<pcl::PointXYZ>::Ptr& curvacloud, int k_neborradius, vector<double>& curvature_thereshold)
{//Mean_curvature.resize(cloud->size(), 1);//初始化矩阵 cloud->size * 1；//Gauss_curvature.resize(cloud->size(), 1);Matrix<double, 6, 6>Coefficient_matrix; //求解方程组系数矩阵Matrix<double, 6, 6>Inverse_Coefficient_matrix;		//系数矩阵的逆矩阵Matrix<double, 6, 1>Value_matrix;					//方程组右值矩阵pcl::KdTreeFLANN<pcl::PointXYZ> kdtree;kdtree.setInputCloud(cloud);vector<int> pointIdRadiusSearch;             //点号vector<float> pointRadiusSquaredDistance;	 //距离for (int i = 0; i < cloud->points.size(); i++){kdtree.nearestKSearch(cloud->points[i], k_neborradius, pointIdRadiusSearch, pointRadiusSquaredDistance);double Target_x = 0, Target_y = 0, Target_z = 0;//数据点的坐标值Coefficient_matrix.setZero();					//初始化矩阵元素全为0Inverse_Coefficient_matrix.setZero();Value_matrix.setZero();double a = 0, b = 0, c = 0, d = 0, e = 0, f = 0;//二次曲面方程系数double u = 0, v = 0;							//二次曲面参数方程参数double E = 0, G = 0, F = 0, L = 0, M = 0, N = 0;//曲面第一、第二基本量Target_x = cloud->points[i].x;Target_y = cloud->points[i].y;Target_z = cloud->points[i].z;double Meancurvature = 0, Gausscurvature = 0;	//平均曲率、高斯曲率//遍历数据点i的邻域点for (int j = 0; j < pointIdRadiusSearch.size(); j++){Coefficient_matrix(0, 0) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x);Coefficient_matrix(0, 1) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(0, 2) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(0, 3) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x);Coefficient_matrix(0, 4) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(0, 5) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x);Coefficient_matrix(1, 1) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(1, 2) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(1, 3) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(1, 4) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(1, 5) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(2, 2) += (cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(2, 3) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(2, 4) += (cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(2, 5) += (cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(3, 3) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x);Coefficient_matrix(3, 4) += (cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(3, 5) += (cloud->points[pointIdRadiusSearch[j]].x);Coefficient_matrix(4, 4) += (cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(4, 5) += (cloud->points[pointIdRadiusSearch[j]].y);Coefficient_matrix(5, 5) += 1;Value_matrix(0, 0) += (cloud->points[pointIdRadiusSearch[j]].z * cloud->points[pointIdRadiusSearch[j]].x * cloud->points[pointIdRadiusSearch[j]].x);Value_matrix(1, 0) += (cloud->points[pointIdRadiusSearch[j]].z * cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].x);Value_matrix(2, 0) += (cloud->points[pointIdRadiusSearch[j]].z * cloud->points[pointIdRadiusSearch[j]].y * cloud->points[pointIdRadiusSearch[j]].y);Value_matrix(3, 0) += (cloud->points[pointIdRadiusSearch[j]].z * cloud->points[pointIdRadiusSearch[j]].x);Value_matrix(4, 0) += (cloud->points[pointIdRadiusSearch[j]].z * cloud->points[pointIdRadiusSearch[j]].y);Value_matrix(5, 0) += (cloud->points[pointIdRadiusSearch[j]].z);}Coefficient_matrix(0, 0) += Target_x * Target_x * Target_x * Target_x;Coefficient_matrix(0, 1) += Target_x * Target_x * Target_x * Target_y;Coefficient_matrix(0, 2) += Target_x * Target_x * Target_y * Target_y;Coefficient_matrix(0, 3) += Target_x * Target_x * Target_x;Coefficient_matrix(0, 4) += Target_x * Target_x * Target_y;Coefficient_matrix(0, 5) += Target_x * Target_x;Coefficient_matrix(1, 1) += Target_x * Target_x * Target_y * Target_y;Coefficient_matrix(1, 2) += Target_x * Target_y * Target_y * Target_y;Coefficient_matrix(1, 3) += Target_x * Target_x * Target_y;Coefficient_matrix(1, 4) += Target_x * Target_y * Target_y;Coefficient_matrix(1, 5) += Target_x * Target_y;Coefficient_matrix(2, 2) += Target_y * Target_y * Target_y * Target_y;Coefficient_matrix(2, 3) += Target_x * Target_y * Target_y;Coefficient_matrix(2, 4) += Target_y * Target_y * Target_y;Coefficient_matrix(2, 5) += Target_y * Target_y;Coefficient_matrix(3, 3) += Target_x * Target_x;Coefficient_matrix(3, 4) += Target_x * Target_y;Coefficient_matrix(3, 5) += Target_x;Coefficient_matrix(4, 4) += Target_y * Target_y;Coefficient_matrix(4, 5) += Target_y;Coefficient_matrix(5, 5) += 1;Value_matrix(0, 0) += Target_x * Target_x * Target_z;Value_matrix(1, 0) += Target_z * Target_x * Target_y;Value_matrix(2, 0) += Target_z * Target_y * Target_y;Value_matrix(3, 0) += Target_z * Target_x;Value_matrix(4, 0) += Target_z * Target_y;Value_matrix(5, 0) += Target_z;Coefficient_matrix(1, 0) = Coefficient_matrix(0, 1);Coefficient_matrix(2, 0) = Coefficient_matrix(0, 2);Coefficient_matrix(2, 1) = Coefficient_matrix(1, 2);Coefficient_matrix(3, 0) = Coefficient_matrix(0, 3);Coefficient_matrix(3, 1) = Coefficient_matrix(1, 3);Coefficient_matrix(3, 2) = Coefficient_matrix(2, 3);Coefficient_matrix(4, 0) = Coefficient_matrix(0, 4);Coefficient_matrix(4, 1) = Coefficient_matrix(1, 4);Coefficient_matrix(4, 2) = Coefficient_matrix(2, 4);Coefficient_matrix(4, 3) = Coefficient_matrix(3, 4);Coefficient_matrix(5, 0) = Coefficient_matrix(0, 5);Coefficient_matrix(5, 1) = Coefficient_matrix(1, 5);Coefficient_matrix(5, 2) = Coefficient_matrix(2, 5);Coefficient_matrix(5, 3) = Coefficient_matrix(3, 5);Coefficient_matrix(5, 4) = Coefficient_matrix(4, 5);Inverse_Coefficient_matrix = Coefficient_matrix.inverse();for (int l = 0; l < 6; ++l){a += Inverse_Coefficient_matrix(0, l) * Value_matrix(l, 0);}for (int l = 0; l < 6; ++l){b += Inverse_Coefficient_matrix(1, l) * Value_matrix(l, 0);}for (int l = 0; l < 6; ++l){c += Inverse_Coefficient_matrix(2, l) * Value_matrix(l, 0);}for (int l = 0; l < 6; ++l){d += Inverse_Coefficient_matrix(3, l) * Value_matrix(l, 0);}for (int l = 0; l < 6; ++l){e += Inverse_Coefficient_matrix(4, l) * Value_matrix(l, 0);}for (int l = 0; l < 6; ++l){f += Inverse_Coefficient_matrix(5, l) * Value_matrix(l, 0);}//根据所求曲面方程的系数计算曲面第一第二基本量u = 2 * a * cloud->points[i].x + b * cloud->points[i].y + d;v = 2 * c * cloud->points[i].y + b * cloud->points[i].x + e;E = 1 + u * u;F = u * v;G = 1 + v * v;double u_v = sqrt(1 + u * u + v * v);L = (2 * a) / u_v;M = b / u_v;N = (2 * c) / u_v;//Gauss_curvatureGausscurvature = (L * N - M * M) / (E * G - F * F);//Mean_curvatureMeancurvature = (E * N - 2 * F * M + G * L) / (2 * E * G - 2 * F * F);curvature_thereshold.push_back(abs(Meancurvature));cout << "The " << i << " point, Gauss_curvature = " << Gausscurvature << "Mean_curvature = " << Meancurvature << endl;}
}int main(int argc, char** argv)
{pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);pcl::PointCloud<pcl::PointXYZ>::Ptr cloudout(new pcl::PointCloud<pcl::PointXYZ>);pcl::io::loadPCDFile<pcl::PointXYZ>("bunny.pcd", *cloud);vector<double> Meancurvature;surfacefitting(cloud, cloudout, 10, Meancurvature);vector<double> curvatureTem = Meancurvature;sort(curvatureTem.begin(), curvatureTem.end());int radio = Meancurvature.size() * 0.8;  //控制提取比列cout << radio;for (int i = 0; i < cloud->size(); i++){if (Meancurvature[i] > curvatureTem[radio]){pcl::PointXYZ pt;pt.x = cloud->points[i].x;pt.y = cloud->points[i].y;pt.z = cloud->points[i].z;cloudout->push_back(pt);}}pcl::io::savePCDFile("extract.pcd", *cloudout);boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer(new pcl::visualization::PCLVisualizer("Point Cloud visual"));int showBefore = 0;int showAfter = 1;viewer->createViewPort(0, 0, 0.5, 1, showBefore);viewer->createViewPort(0.5, 0, 1, 1, showAfter);viewer->setBackgroundColor(0, 0, 0, showBefore);viewer->setBackgroundColor(0.05, 0, 0, showAfter);//原始点云pcl::visualization::PointCloudColorHandlerCustom<pcl::PointXYZ> src_h(cloud, 0, 255, 0);//输出点云pcl::visualization::PointCloudColorHandlerCustom<pcl::PointXYZ> after_sac(cloudout, 0, 0, 255);viewer->setBackgroundColor(0, 0, 0);viewer->addPointCloud(cloud, src_h, "source cloud", showBefore);viewer->addPointCloud(cloudout, after_sac, "target cloud1", showAfter);while (!viewer->wasStopped()){viewer->spinOnce(100);boost::this_thread::sleep(boost::posix_time::microseconds(10000));}}

四、最终结果

这篇关于二次曲面拟合计算点云法曲率、主曲率以及平均曲率（手写版）的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/701323。

相关文章

计算绕原点旋转某角度后的点的坐标

计算绕原点旋转某角度后的点的坐标

问题： A点（x, y）按顺时针旋转 theta 角度后点的坐标为A1点(x1,y1)  ，求x1 y1坐标用（x，y）和 theta 来表示 方法一： 设 OA 向量和x轴的角度为 alpha ， 那么顺时针转过 theta后 ，OA1 向量和x轴的角度为 (alpha - theta) 。 使用圆的参数方程来表示点坐标。A的坐标可以表示为： \[\left\{ {\begin{ar
阅读更多...
19.手写Spring AOP

19.手写Spring AOP

1.Spring AOP顶层设计 2.Spring AOP执行流程 下面是代码实现 3.在 application.properties中增加如下自定义配置: #托管的类扫描包路径#scanPackage=com.gupaoedu.vip.demotemplateRoot=layouts#切面表达式expression#pointCut=public .* com.gupaoedu
阅读更多...
17.用300行代码手写初体验Spring V1.0版本

17.用300行代码手写初体验Spring V1.0版本

1.1.课程目标 1、了解看源码最有效的方式，先猜测后验证，不要一开始就去调试代码。 2、浓缩就是精华，用 300行最简洁的代码 提炼Spring的基本设计思想。 3、掌握Spring框架的基本脉络。 1.2.内容定位 1、 具有1年以上的SpringMVC使用经验。 2、 希望深入了解Spring源码的人群，对 Spring有一个整体的宏观感受。 3、 全程手写实现SpringM
阅读更多...
【云计算 复习】第1节 云计算概述和 GFS + chunk

【云计算 复习】第1节 云计算概述和 GFS + chunk

一、云计算概述 1.云计算的商业模式 （1）软件即服务（SaaS） 有些景区给游客提供烧烤场地，游客需要自己挖坑或者砌烧烤台，然后买肉、串串、烧烤。 （2）平台即服务（PaaS） 有些景区给游客提供烧烤场地，同时搭建好烧烤台，游客只需要自己带食材和调料、串串、烧烤。 （3）基础设施即服务（IaaS） 有些景区给游客提供烧烤场地，同时搭建好烧烤台，还有专门的厨师来烧烤，用户不需要关心前面的所有
阅读更多...
神经网络第四篇：推理处理之手写数字识别

神经网络第四篇：推理处理之手写数字识别

到目前为止，我们已经介绍完了神经网络的基本结构，现在用一个图像识别示例对前面的知识作整体的总结。本专题知识点如下： MNIST数据集图像数据转图像神经网络的推理处理批处理  MNIST数据集          mnist数据图像 MNIST数据集由0到9的数字图像构成。像素取值在0到255之间。每个图像数据都相应地标有“7”、“2”、“1”等数字标签。MNIST数据集中，
阅读更多...
什么是dB？dBm、dBc、dBi、dBd怎么计算，有什么区别？

什么是dB？dBm、dBc、dBi、dBd怎么计算，有什么区别？

什么是dB？dBm、dBc、dBi、dBd怎么计算，有什么区别？ 引言 在电子工程、通信和音频领域，dB（分贝）是一个常见的术语。许多人刚接触时可能会感到困惑，因为它不仅仅是一个简单的单位，还有多种不同的形式，如dBm、dBc、dBi和dBd。这篇文章将详细解释这些概念，并介绍如何计算它们，帮助初学者更好地理解和应用。 什么是dB？ dB，即分贝，是一种表示两个数值比值的对数单位。分贝的基
阅读更多...
基于NURBS曲线的数据拟合算法matlab仿真

基于NURBS曲线的数据拟合算法matlab仿真

目录 1.程序功能描述 2.测试软件版本以及运行结果展示 3.核心程序 4.本算法原理 4.1NURBS曲线基础 4.2 数据拟合原理 5.完整程序 1.程序功能描述        基于NURBS曲线的数据拟合算法,非均匀有理B样条（Non-Uniform Rational B-Splines，简称NURBS）曲线是一种强大的数学工具，广泛应用于计算机图形学、CAD/CA
阅读更多...
php字符串计算汉字、中英文数字个数

php字符串计算汉字、中英文数字个数

$str = '123abcDEF测试的事发地点';$length = strlen(preg_replace('/[\x00-\x7F]/', '', $str));$arr['en'] = strlen( $str) - $length; //（非中文）$arr['cn'] = intval($length / 3); // 编码GBK，除以2 （中文）print_r($
阅读更多...
计算广告：第四章——合约广告

计算广告：第四章——合约广告

计算广告：第四章——合约广告 一、广告位合约 二、受众定向 1、受众定向方法概览 2、 受众定向标签体系 三、展示量合约 1、流量预测 2、流量塑性 3、在线分配 包括按 CPM 计费的展示量合约广告和按 CPT 结算的广告位合约。   一、广告位合约 按CPT结算广告位合约 缺点：无法做到按受众类型投放广告，无法进行深入的优化效果 优点：强曝光属性带来品牌冲击，或
阅读更多...
计算广告：第三章——在线广告产品概览

计算广告：第三章——在线广告产品概览

第三章——在线广告产品概览 一、商业产品的设计原则 二、需求方层级组织及接口 二、供给方管理接口 （1）合约广告产品——主要服务于后续效果不宜直接衡量的品牌类广告主 按时段售卖的CPT广告按约定展示量售卖的CPM广告   （2）竞价广告产品 其形式主要是搜索广告，其产品形式为对搜索关键词的竞价。这种广告拓展到站外广告时，演变为了对页面关键词或者用户标签竞价的产品形式，也就是
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2