精华内容
下载资源
问答
  • 2022-04-12 13:25:49

    点云法线估计原理参考:
    [PCL]2 点云法向量计算NormalEstimation
    点云法向量估计原理及应用PCL
    点云曲率计算原理参考:
    二次曲面拟合求点云的曲率
    (曲率系列3:)PCL:PCL库中的两种曲率表示方法pcl::NormalEstimation和PrincipalCurvaturesEstimation
    代码实现:

    /**
     * @description:	最小二乘法拟合平面
     * @param cloud		输入点云
     * @return			平面方程系数
     */
    std::vector<double> FitPlaneByLeastSquares(pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud)
    {
    	if (cloud->points.size() < 3)	return {};
    
    	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 < cloud->points.size(); i++)
    	{
    		meanX += cloud->points[i].x;
    		meanY += cloud->points[i].y;
    		meanZ += cloud->points[i].z;
    
    		meanXX += cloud->points[i].x * cloud->points[i].x;
    		meanYY += cloud->points[i].y * cloud->points[i].y;
    		meanZZ += cloud->points[i].z * cloud->points[i].z;
    
    		meanXY += cloud->points[i].x * cloud->points[i].y;
    		meanXZ += cloud->points[i].x * cloud->points[i].z;
    		meanYZ += cloud->points[i].y * cloud->points[i].z;
    	}
    	meanX /= cloud->points.size();
    	meanY /= cloud->points.size();
    	meanZ /= cloud->points.size();
    	meanXX /= cloud->points.size();
    	meanYY /= cloud->points.size();
    	meanZZ /= cloud->points.size();
    	meanXY /= cloud->points.size();
    	meanXZ /= cloud->points.size();
    	meanYZ /= cloud->points.size();
    
    	Eigen::Matrix3d m;
    	m(0, 0) = meanXX - meanX * meanX; m(0, 1) = meanXY - meanX * meanY; m(0, 2) = meanXZ - meanX * meanZ;
    	m(1, 0) = meanXY - meanX * meanY; m(1, 1) = meanYY - meanY * meanY; m(1, 2) = meanYZ - meanY * meanZ;
    	m(2, 0) = meanXZ - meanX * meanZ; m(2, 1) = meanYZ - meanY * meanZ; m(2, 2) = meanZZ - meanZ * meanZ;
    	Eigen::EigenSolver<Eigen::Matrix3d> PlMat(m * cloud->points.size());
    	Eigen::Matrix3d eigenvalue = PlMat.pseudoEigenvalueMatrix();
    	Eigen::Matrix3d eigenvector = PlMat.pseudoEigenvectors();
    
    	double v1 = eigenvalue(0, 0), v2 = eigenvalue(1, 1), v3 = eigenvalue(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 = eigenvector(0, minNumber), b = eigenvector(1, minNumber), c = eigenvector(2, minNumber), d = -(a * meanX + b * meanY + c * meanZ);
    
    	return{ a / sqrt(a*a + b*b + c*c), b / sqrt(a*a + b*b + c*c) , c / sqrt(a*a + b*b + c*c), std::min(v1,std::min(v2,v3)) / (v1 + v2 + v3) };
    }
    
    /**
     * @description:	最小二乘法拟合二次曲面
     * @param cloud		输入点云
     * @param point		给定点
     * @return			二次曲面给定点初的高斯曲率
     */
    double FitQuadricByLeastSquares(pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, pcl::PointXYZ& point)
    {
    	if (cloud->points.size() < 6)	return{};
    
    	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;	//二次曲面基本量
    
    	double mean_curvature = 0, guass_curvature = 0;		//平均曲率和高斯曲率
    
    	Eigen::Matrix<double, 6, 6> Q; //线性方程组系数矩阵
    	Eigen::Matrix<double, 6, 1> B; //线性方程组右值矩阵
    
    	Q.setZero();
    	B.setZero();
    
    	for (int i = 0; i < cloud->points.size(); i++)
    	{
    		Q(0, 0) += pow(cloud->points[i].x, 4);
    		Q(1, 0) = Q(0, 1) += pow(cloud->points[i].x, 3)*cloud->points[i].y;
    		Q(2, 0) = Q(0, 2) += pow(cloud->points[i].x*cloud->points[i].y, 2);
    		Q(3, 0) = Q(0, 3) += pow(cloud->points[i].x, 3);
    		Q(4, 0) = Q(0, 4) += pow(cloud->points[i].x, 2)*cloud->points[i].y;
    		Q(5, 0) = Q(0, 5) += pow(cloud->points[i].x, 2);
    		Q(1, 1) += pow(cloud->points[i].x*cloud->points[i].y, 2);
    		Q(2, 1) = Q(1, 2) += pow(cloud->points[i].y, 3)*cloud->points[i].x;
    		Q(3, 1) = Q(1, 3) += pow(cloud->points[i].x, 2)*cloud->points[i].y;
    		Q(4, 1) = Q(1, 4) += pow(cloud->points[i].y, 2)*cloud->points[i].x;
    		Q(5, 1) = Q(1, 5) += cloud->points[i].x*cloud->points[i].y;
    		Q(2, 2) += pow(cloud->points[i].y, 4);
    		Q(3, 2) = Q(2, 3) += pow(cloud->points[i].y, 2)*cloud->points[i].x;
    		Q(4, 2) = Q(2, 4) += pow(cloud->points[i].y, 3);
    		Q(5, 2) = Q(2, 5) += pow(cloud->points[i].y, 2);
    		Q(3, 3) += pow(cloud->points[i].x, 2);
    		Q(4, 3) = Q(3, 4) += cloud->points[i].x*cloud->points[i].y;
    		Q(5, 3) = Q(3, 5) += cloud->points[i].x;
    		Q(4, 4) += pow(cloud->points[i].y, 2);
    		Q(5, 4) = Q(4, 5) += cloud->points[i].y;
    		Q(5, 5) += 1;
    
    		B(0, 0) += pow(cloud->points[i].x, 2)*cloud->points[i].z;
    		B(1, 0) += cloud->points[i].x*cloud->points[i].y*cloud->points[i].z;
    		B(2, 0) += pow(cloud->points[i].y, 2)*cloud->points[i].z;
    		B(3, 0) += cloud->points[i].x*cloud->points[i].z;
    		B(4, 0) += cloud->points[i].y*cloud->points[i].z;
    		B(5, 0) += cloud->points[i].z;
    
    		Eigen::Matrix<double, 6, 6> Q_inverse = Q.inverse();
    		
    		//求解二次曲面方程系数z=ax^2+bxy+cy^2+dx+ey+f
    		for (size_t j = 0; j < 6; ++j)
    		{
    			a += Q_inverse(0, j)*B(j, 0);
    			b += Q_inverse(1, j)*B(j, 0);
    			c += Q_inverse(2, j)*B(j, 0);
    			d += Q_inverse(3, j)*B(j, 0);
    			e += Q_inverse(4, j)*B(j, 0);
    			f += Q_inverse(5, j)*B(j, 0);
    		}		
    	}
    
    	//求解二次曲面参数方程系数u=dz/dx,v=dz/dy
    	u = 2 * a * point.x + b * point.y + d;
    	v = 2 * c * point.y + b * point.x + e;
    
    	//求解二次曲面第一基本量 E = u^2, F = uv, G = v^2
    	E = 1 + u * u;
    	F = u * v;
    	G = 1 + v * v;
    
    	//求解二次曲面第二基本量 L = p*n, M = r*n, N = q*n; p = d^2z/dx^2, q = d^2z/dy^2, r = d^2z/dxdy;n = (u×v) / |u×v|
    	L = 2 * a / sqrt(1 + u * u + v * v);
    	M = b / sqrt(1 + u * u + v * v);
    	N = 2 * c / sqrt(1 + u * u + v * v);
    
    	// 高斯曲率
    	guass_curvature = (L * N - M * M) / (E * G - F * F);
    
    	// 平均曲率
    	mean_curvature = (E * N - 2 * F * M + G * L) / (2 * E * G - 2 * F * F);
    
    	return guass_curvature;
    }
    
    /**
     * @description:	法线估计
     * @param cloud		输入点云
     * @param normals	法向量和曲率
     * @param nr_k		k邻近点数
     */
    void normalestimation(pcl::PointCloud<pcl::PointXYZ>::Ptr& cloud, pcl::PointCloud <pcl::Normal>::Ptr normals, int nr_k)
    {
    	pcl::KdTreeFLANN<pcl::PointXYZ> kdtree;
    	kdtree.setInputCloud(cloud);
    	std::vector<int> pointsIdx(nr_k);          
    	std::vector<float> pointsDistance(nr_k);   
    
    	for (size_t i = 0; i < cloud->size(); ++i)
    	{
    		kdtree.nearestKSearch(cloud->points[i], nr_k, pointsIdx, pointsDistance);    //k近邻搜索
    		pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_temp(new pcl::PointCloud<pcl::PointXYZ>);
    		pcl::copyPointCloud(*cloud, pointsIdx, *cloud_temp);
    	
    		std::vector<double> norm = FitPlaneByLeastSquares(cloud_temp);	//计算法向量
    		//double curvature = FitQuadricByLeastSquares(cloud_temp, cloud->points[i]);		//计算曲率
    
    		pcl::Normal p(norm[0], norm[1], norm[2]);
    		p.curvature = norm[3];
    		normals->push_back(p);
    	}
    }
    

    代码传送门:https://github.com/taifyang/PCL-algorithm

    更多相关内容
  • 点云曲率计算

    2021-03-23 01:30:18
    点云数据计算曲率,使用C++编写的程序,计算每个点的曲率 点云数据计算曲率,使用C++编写的程序,计算每个点的曲率
  • 离散点曲率计算离散点曲率计算离散点曲率计算离散点曲率计算
  • [K,H,P1,P2] = surfature(X,Y,Z) 返回曲面的高斯曲率 (K)、平均曲率 (H) 和主曲率 (P1,P2)。 输入 (X,Y,Z) 是对应于被分析表面的二维数组。 例子[X,Y,Z] = 峰值; [K,H,P1,P2] = 表面(X,Y,Z); 冲浪(X,Y,Z,H,'...
  • NURBS曲线曲率计算

    2018-04-10 10:38:00
    计算NURBS曲线在给定的节点处的曲率,需要配合NURBS工具箱使用。
  • matlab开发-三角曲面3的曲率计算。在三维三角形网格上计算主曲率。
  • 用于计算三角形网格上主曲率的函数。 曲率近似值基于局部 (N=1) 邻域元素和顶点。 注意:曲率方向未正确计算。 一旦实施,将发布更新版本。 请注意,相邻三角形很少的顶点处的计算,因此相邻顶点很少,被扩展到更...
  • 计算nurbs曲率NURBS曲线的曲率计算
  • 空间曲线的曲率计算方法(附代码)

    千次阅读 2022-02-18 21:25:56
    空间曲线的曲率计算方法,附代码

    一、空间曲线的曲率计算方法

      方法一,参考张学东的论文《空间曲线的曲率计算方法》:
    在这里插入图片描述
      方法二,参考博文:空间曲线的弧长与曲率
    在这里插入图片描述

    二、空间曲线的曲率计算matlab代码

    %{
    Function: calculate_curvature_of_spatial_curve
    Description: 计算空间曲线的曲率
    Input: 速度向量v,加速度向量a
    Output: 空间曲线的曲率k,求解结果状态sta(sta = 0表示求解失败,sta = 1表示求解成功)
    Author: Marc Pony(marc_pony@163.com)
    %}
    function [k, sta] = calculate_curvature_of_spatial_curve(v, a)
    
    normV = norm(v, 2);
    normA = norm(a, 2);
    dotVA = dot(v, a);
    temp = (normA * normV - dotVA) * (normA * normV + dotVA);
    if temp >= 0.0 && normV > eps
        k = sqrt(temp) / normV / normV / normV;
        sta = 1;
    else
        k = 0.0;
        sta = 0;
    end
    end
    
    %{
    Function: calculate_curvature_of_spatial_curve2
    Description: 计算空间曲线的曲率
    Input: 速度向量v,加速度向量a
    Output: 空间曲线的曲率k,求解结果状态sta(sta = 0表示求解失败,sta = 1表示求解成功)
    Author: Marc Pony(marc_pony@163.com)
    %}
    function [k, sta] = calculate_curvature_of_spatial_curve2(v, a)
    
    normV = norm(v, 2);
    if normV > eps
        k = norm(cross(v, a), 2) / normV / normV / normV;
        sta = 1;
    else
        k = 0.0;
        sta = 0;
    end
    
    end
    
    clc
    clear
    close all
    
    %三次Bezier曲线C(u) = (1-u)^3 * P0 + 3*u*(1-u)^2*P1 + 3*u^2*(1-u)*P2 + u^3*P3
    %                  = (3*P1 - P0 - 3*P2 + P3)*u^3 + (3*P0 - 6*P1 + 3*P2)*u^2 + (3*P1 - 3*P0)*u + P0
    
    syms P0 P1 P2 P3 u
    collect(expand((1-u)^3 * P0 + 3*u*(1-u)^2*P1 + 3*u^2*(1-u)*P2 + u^3*P3 ))
    
    P0 = [0, 0, 0];
    P1 = [10, 30, 0];
    P2 = [20, 5, 0];
    P3 = [30, 30, 0];
    
    n = 300;
    u = linspace(0, 1, n);
    k = zeros(n, 1);
    k2 = zeros(n, 1);
    sta = zeros(n, 1);
    sta2 = zeros(n, 1);
    p = zeros(n, 3);
    v = zeros(n, 3);
    a = zeros(n, 3);
    for i = 1 : n
        p(i, :) = (3*P1 - P0 - 3*P2 + P3)*u(i)^3 + (3*P0 - 6*P1 + 3*P2)*u(i)^2 + (3*P1 - 3*P0)*u(i) + P0;
        v(i, :) = 3 * (3*P1 - P0 - 3*P2 + P3)*u(i)^2 + 2*(3*P0 - 6*P1 + 3*P2)*u(i) + (3*P1 - 3*P0);
        a(i, :) = 6 * (3*P1 - P0 - 3*P2 + P3)*u(i) + 2*(3*P0 - 6*P1 + 3*P2);
        [k(i), sta(i)] = calculate_curvature_of_spatial_curve(v(i, :), a(i, :));
        [k2(i), sta2(i)] = calculate_curvature_of_spatial_curve2(v(i, :), a(i, :));
    end
    error = sum(abs(k-k2))
    
    R = 1./ k;  %曲率半径
    
    figure
    subplot(3, 1, 1)
    plot3(p(:, 1), p(:, 2), p(:, 3), '-')
    hold on
    plot3([P0(1), P1(1), P2(1), P3(1)],[P0(2), P1(2), P2(2), P3(2)],[P0(3), P1(3), P2(3), P3(3)], 'o-')
    ylabel('Bezier曲线')
    view([0 0 1])
    subplot(3, 1, 2)
    plot(u, k, '-')
    ylabel('曲率')
    subplot(3, 1, 3)
    plot(u, R, '-')
    ylabel('曲率半径')
    

    在这里插入图片描述

    展开全文
  • 图像工程中,离散的点云模型表面曲率的精确计算非常重要。笔者在Levin 的MLS(Moving Least-Square)表面的基础上,将原始三维点云模型的表面投影到一个MLS表面,然后直接从MLS表面计算点云曲率,并将该方法应用在...
  • 对于做飞思卡尔比赛的同学来说,智能车在赛道上的赛道曲率是相当重要,对于智能车车的控制起到至关重要的参考
  • 曲率属性近年来在构造识别和解释上得到了迅速的发展和应用,而曲率计算通常使用3×3网格单元对局部曲面作最小二乘法逼近,由Roberts推导并给出了具体的计算公式,该方法在求取曲率值时带来了大量噪声。采用5×5网格单元...
  • 基于 Bragg光纤光栅的中心波长随应变发生偏移的特点,提出利用三芯光纤的 Bragg光栅进行三维曲率计算的理论.分析了 3个光纤光栅的中心波长随曲率半径、中性面到中心的距离、偏转角之间的变化关系,并且用 MATLAB进行...
  • 轮廓曲率计算(附python代码)

    千次阅读 2021-11-27 20:03:49
    最近在研究角点的检测,需要用到曲率计算,就研究了一下目前常用的曲率计算,我这里用的是点到切线距离累加和曲率估计方法(RTPDA) 参考了一下重庆大学的一篇博士论文《几种轮廓曲率估计角点检测算法研究》 先附上...

    最近在研究角点的检测,需要用到曲率计算,就研究了一下目前常用的曲率计算,我这里用的是点到切线距离累加和曲率估计方法(RTPDA)
    参考了一下重庆大学的一篇博士论文《几种轮廓曲率估计角点检测算法研究
    先附上链接曲率计算代码

    看起来公式有点多,其实理解起来还是很容易,请大家仔细看下去。

    另外关于最小二乘法,可以参考链接最小二乘法的本质

    还有一些简单的计算曲率的方案,比如“直接利用曲率的计算公式”,或者是“通过三角形外接圆半径求曲率”,这些对于离散点来说,还是有点影响的
    曲率计算

    展开全文
  • 曲率计算过程: 1. 找出所有的公共边,以及包含他们的两个三角形面片。 保存公共边起始-终点顶点编号(VerticeID, VerticeID)、三角形面片(TriaID,TriaID)编号、单位化的公共边向量edgeVector=[vx, vy,v

    参考:

    https://blog.csdn.net/weixin_44210987/article/details/113279727

    https://blog.csdn.net/qq_36686437/article/details/116369956

    https://blog.csdn.net/qq_36686437/article/details/105559280

    https://blog.csdn.net/ModestBean/article/details/89438082

    三维空间中的曲率:三维曲面偏离平面的程度           

                                 

    曲面曲率:
    在曲面上取一点P,曲面在P点的法线为n,过n可以有无限多个剖切平面,每个剖切平面与曲面相交,交线为一条平面曲线。平面结论:圆上弯曲程度相同,任意一点曲率相等,越弯曲曲率越大,直线曲率为0。

    不同的剖切平面上的平面曲线在P点的曲率半径一般是不相等的。

    主曲率:曲面上有无数个不同方向的曲线,曲面上的点不同方向具有不同曲率,其中最大值和最小值为称为主曲率 k1 和k2,极值方向称为主方向。数学上可证名k1和k2互相垂直。

    高斯曲率:两主曲率乘积,反映曲面在不同方向弯曲程度是否相同。高斯曲率为正,为球面。高斯曲率为负双曲面。

    平均曲率:两主曲率算数平均数(k1+k2)/2,反映曲面凹凸程度。平均曲率为正,局部凹。平均曲率为负,局部凸。
     

    曲率计算过程:

    1. 找出所有的公共边,以及包含他们的两个三角形面片。

    保存公共边起始-终点顶点编号(VerticeID, VerticeID)、三角形面片(TriaID,TriaID)编号、单位化的公共边向量edgeVector=[vx, vy,vz] 以及公共边向量长度distances。

    2. 求相邻两个三角形面片法向量的cos夹角beta,并符号化。

    符号化:

    相邻三角形法向量叉乘之后的向量cp,与edgeVector 进行点乘,若结果大于0,sign=1;等于0,sign=0;小于0,sign=-1.

    ### 右手系法则下,

    两个三角形为凸,则cp与edgeVector夹角为0, sign=1;

    两个三角形为凹,则cp与edgeVector夹角为180, sign=-1;
    两个三角形呈一条线,它们的法向量平行,cp==0, sign=0;

    3. 构建曲率公式:

    T =f(edgeVector, beta,sign,distance)

    4. 矩阵分解,求特征向量和特征值。

    特征值最小的为最小曲率, 

    特征值最大的为法向量, 

    特征值第二大的的为最大曲率

                Cmean = (Cmin+Cmax)/2
                Cgauss = Cmin*Cmax

       

    2. input:
            vertices: [nx3]
            faces:[n,3]
            normals:[n,3]
        
        calculate:
            edgeVector:相邻三角形的公共边单位向量。[3,ne], ne为公共边数量
            beta:相邻两个三角形法向量的cos角。[1,ne]
            Tv:曲率公式
        return:
            Umin, Umax, Cmin, Cmax, Cmean, Cgauss
            依次为最小最大切向量,最小最大曲率,平均曲率,高斯曲率。

    import numpy as np
    from numpy import linalg as LA
    import sys
    
    
    '''
    计算平均曲率和高斯曲率
    
    '''
    
    def getCommonEdges(faces):
        '''
        input:
            faces:[nx3]
        return:
            commonEdges: indexes of vertices in common edges, [nx2]
            facePairs:  triangle pairs where they locate in ,[nx2]
        '''
    
        # faces = np.array([[1,2,3],[3,2,4],[5,2,1]])
        faces = faces-1### python 索引 从 0 开始
        numFace = faces.shape[0]
        edgeStart = faces.flatten() 
        edgeEnd = faces[:,[1,2,0]].flatten() 
        edges = np.vstack((edgeStart, edgeEnd)) #[2,n]
        faceId = np.tile(np.arange(numFace),(3,1)).transpose().flatten().tolist()
    
        commonEdges = []
        facePires = []
        numEdges= edges.shape[1]
        for i in range(numEdges):
            curEdge = edges[:,i].tolist()
            # print(curEdge)
            ind1 = np.where(edges[1,:] == curEdge[0])[0]
            ind2 = np.where(edges[0,:] == curEdge[1])[0]
            index = list(set(ind1).intersection(set(ind2)))
            if len(index):
                commonEdges.append(curEdge)
                facePires.append([faceId[i], faceId[index[0]]])
        
        commonEdges=np.array(commonEdges)
        facePires = np.array(facePires)
    
        ###### 去除冗余:edgeStart<edgeEnd ################
        directioned = np.where(commonEdges[:,0]<commonEdges[:,1])[0]
        commonEdges = commonEdges[directioned,:]
        facePires = facePires[directioned,:]
    
        return commonEdges, facePires
    
    
    
    def getCurvature(vertices, faces, normals):
        '''
        input:
            vertices: [nx3]
            faces:[n,3]
            normals:[n,3]
        
        calculate:
            edgeVector:相邻三角形的公共边单位向量。[3,ne], ne为公共边数量
            beta:相邻两个三角形法向量的cos角。[1,ne]
            Tv:曲率公式
        return:
            Umin, Umax, Cmin, Cmax, Cmean, Cgauss
            依次为最小最大切向量,最小最大曲率,平均曲率,高斯曲率。
    
        '''
    
        vertices = np.array(vertices).transpose() #[3,n]
        normals = np.array(normals).transpose() #[3,n]
        numVertices = vertices.shape[1]
    
        ### 1. find common edges and corresponding pairs of triangles
        commonEdges, facePires = getCommonEdges(faces) #[n,2],[n,2]
        ne = commonEdges.shape[0]
        #### normalized edge vector #######
        edgeStart = commonEdges[:,0]
        edgeEnd = commonEdges[:,1]
        edgeVector = vertices[:,edgeEnd] - vertices[:,edgeStart] #[3,ne]
        distances = np.sqrt(np.sum(edgeVector**2, axis = 0)) #[1,ne]
        edgeVector = edgeVector/np.tile(distances,(3,1))
        distances= distances/distances.mean()
    
        ##### 2. cos angle: beta , 具有公共边的两个三角形法向量的夹角 ############
        faceInd1 = facePires[:,0]
        faceInd2 = facePires[:,1]
        dp = np.sum(normals[:,faceInd1] * normals[:,faceInd2], axis=0)
        dp = np.maximum(-1, dp)
        dp = np.minimum(1, dp)
        beta = np.arccos(dp) #[1,ne]
        #### sign: positive or negtive ############
        ### 右手系法则下,两个三角形为凸,则cp与edgeVector夹角为0, sign=1;两个三角形为凹,则cp与edgeVector夹角为180, sign=-1;
        ### 两个三角形一条线,它们的法向量平行,cp==0, sign=1;
        cp = np.cross(normals[:,faceInd1].transpose(), normals[:, faceInd2].transpose()).transpose()#cp: [3,ne]
        sign = np.sign(np.sum((cp * edgeVector), axis=0))
        beta = beta*sign
    
        ###### 3. 构建曲率函数 ##############
        T = np.zeros((3,3,ne))
        for i in range(3):
            for j in range(3):
                T[i,j,:] = np.reshape(edgeVector[i,:]*edgeVector[j,:], (1,1,ne))
                T[j,i,:] = T[i,j,:] 
        ###最后的曲率公式如下 ###
        T = T * np.tile(np.reshape(distances * beta, (1,1,ne)), (3,3,1))
    
    
        ###### 4. 构建关于所有顶点的矩阵(一个顶点一个通道),并将上述结算结果赋给公共边所含顶点 ###################
        Tv = np.zeros((3,3,numVertices))
        w = np.zeros((1,1,numVertices))
        for k in range(ne): 
            Tv[:,:,edgeStart[k]] = Tv[:,:,edgeStart[k]] + T[:,:,k]
            Tv[:,:,edgeEnd[k]] = Tv[:,:,edgeEnd[k]] + T[:,:,k]
            w[:,:,edgeStart[k]] = w[:,:,edgeStart[k]] + 1
            w[:,:,edgeEnd[k]] = w[:,:,edgeEnd[k]] + 1
        # eps = eps(1) = 2.2204e-16;eps为系统运算时计算机允许取到的最小值
        w = np.where(w<sys.float_info.epsilon, 1, w)
        Tv = Tv/np.tile(w, (3,3,1))
    
        # ###### 5. smoothing ##############
        # for x in range(3):
        #     for y in range(3):
        #         a = Tv[x,y,:]
        #         a = meshSmoothing(faces,vertices,a(:),options)
        #         Tv[x,y,:] = a
    
        ###### 6. 矩阵分解:求特征向量eigenvectors and 特征值eigenvalues  ###############
        U = np.zeros((3,3,numVertices))
        D = np.zeros((3,numVertices))
        for k in range(numVertices):
            [d, u] = LA.eig(Tv[:,:,k])
            d = np.real(d)
            ### sort: 按照特征值的绝对值从小到大排序 ######
            sortedInd = sorted(range(len(d)), key=lambda k: abs(d)[k])
            D[:,k] = d[sortedInd]
            U[:,:,k] = np.real(u[:,sortedInd])
        
        Umin =np.squeeze(U[:,2,:]) #### [3,n]
        Umax = np.squeeze(U[:,1,:]) #### [3,n]
        Cmin = D[1,:].transpose() #### [1,n]
        Cmax = D[2,:].transpose() #### [1,n]
        Cmean = (Cmin+Cmax)/2
        Cgauss = Cmin*Cmax
        Normal = np.squeeze(U[:,0,:])
    
        return Umin, Umax, Cmin, Cmax, Cmean, Cgauss
    
    
    
    if __name__=='__main__':
        from readData.objParser import *
        modelFile ='building2-1/components/pillar1.obj'
        
        ###### 数据解析 ########################
        objParser = OBJ_PARSER(modelFile)
        faces = objParser.get_faces()
        faces = np.array(faces)[:,[0,2,4]]
        vertices = objParser.get_vertices()
        vertices = np.array(vertices).astype(np.float64)
        normals = objParser.get_normals()
        ########### 计算曲率 #####################
        Umin, Umax, Cmin, Cmax, Cmean, Cgauss = getCurvature(vertices, faces, normals)
    
    
    
    
    
    
    

    展开全文
  • 曲率第!+卷第(期(""(年&月塔里木农垦大学学报<=>9?5@=AB59CDE?C3F9GC/0=AHI9CJ>...空间曲线的曲率计算方法张学东(塔里木农垦大学农业工程学院,新疆阿拉尔’+))"")通常我们知道如何求一个...
  • Matlab曲率、平均曲率计算

    千次阅读 2021-04-05 17:23:41
    Matlab平均曲率计算计算公式代码使用实列end 计算公式 代码 整理为了一个函数. // 输出的是平均曲率,曲率序列 function [mean1,k2] = Mean_curvature(x0,y0) h = abs(diff([x0(2), x0(1)])); %一阶导 yapp1 = ...
  • 曲率计算公式推导

    千次阅读 2021-01-17 12:40:34
    曲率半径公式推导曲率(k):描述曲线下降长度随角度变化,${\rm{k}} = \mathop {\lim }\limits_{\alpha \to 0} \left| {\frac{{\Delta \alpha }}{{\Delta s}}} \right|$$R = \frac{1}{k} = \frac{{{{\left[ {1 + {{\...
  • PCL:从法线计算到曲率计算并可视化

    千次阅读 多人点赞 2020-04-28 16:53:43
    // 建立主曲率计算 pcl::PrincipalCurvaturesEstimation principalCurvaturesEstimation; // 提供原始点云(没有法线) principalCurvaturesEstimation.setInputCloud (cloud); // 为点云提供法线 ...
  • 简单实用曲率计算工具
  • nurbs曲线的曲率进行计算,注意: ①曲率半径为曲率的倒数 ②如果是离散点,先用polyfit和polyval拟合出曲线
  • 离散点的曲率计算

    千次阅读 2020-08-08 19:16:26
    针对曲线上某个点的切线方向角对弧长的转动率,而曲率的倒数就是曲率半径 K  =  ∣Δθl∣  =  ∣1r∣  K\; =\; \left| \frac{\Delta \theta }{l} \right|\; =\; \left| \frac{1}{r} \right|\; K=∣∣∣∣...
  • 本发明涉及工程零件表面点云测量领域,更具体的,涉及一种点云主曲率计算方法。背景技术:三维数据测量现已广泛应用于复杂型面零件的加工质量检测、运行状态监测及快速检修过程中,以及时进行外形评价并做出决策,...
  • matlab算法计算三维散乱点云的曲率,包括主曲率,高斯曲率和平均曲率 matlab算法计算三维散乱点云的曲率,包括主曲率,高斯曲率和平均曲率
  • 曲率计算公式

    2017-03-22 09:59:25
    曲率计算公式以及推导过程
  • 之前在学习法向量的计算时曾经说过,PCA算法最初其实是一种评估曲率的方法,表示一个点的曲率有很多种,但是一般他们都需要使用三角网格来进行表示,文献[1]中的作者通过定义高斯曲率KKK以及平均曲率HHH来表示一个...
  • I know the edge detection problem has been posted before (in Java: Count the number of objects in an Image, language independent: Image edge detection), but I want to know how to implement it in pytho...
  • 用于飞思卡尔智能车比赛摄像头组的曲率计算

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 12,615
精华内容 5,046
关键字:

曲率计算

友情链接: rectangular.zip