三角网格的顶点曲率计算（平均曲率和高斯曲率）

参考：

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, facePiresdef 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 trianglescommonEdges, 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]] + 1w[:,:,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)/2Cgauss = Cmin*CmaxNormal = np.squeeze(U[:,0,:])return Umin, Umax, Cmin, Cmax, Cmean, Cgaussif __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)

三角网格的顶点曲率计算（平均曲率和高斯曲率）相关推荐

三维网格上的曲率计算
[1] 基本概念曲率是用来反映几何体的弯曲程度,定性的看,弯曲的越厉害,该部分的曲率越大. 平均曲率.主曲率和高斯曲率是曲率的三个基本要素. 平均曲率:是空间上曲面上某一点任意两个相互垂直的正交曲率 ...
数字几何处理作业1：编程实现三角网格上高斯曲率和平均曲率的计算编程部分
三.编程 1.代码用的是中国科大傅孝明老师的框架:框架下载及配置运行 (1)在哪儿添加代码梳理框架的结构后,在MeshViewerWidget.中添加求解曲率的函数,并在MainViewerWid ...
三角网格模型及基于RBF隐曲面方程求解的曲面重建
资料来源:径向基函数和神经网络技术在逆向工程中的应用研究(博士论文:王宏涛) RBF神经网络模型 RBF神经网络起源于数值分析中多变量插值的RBF方法,1988年Broomhead等人首先将该算法应用 ...
Open3D 点云曲率计算
文章目录一.简介二.代码实现三.相关函数(注意) 四.实现效果参考资料一.简介由于在Open3D中也没发现相关函数,就自己动手丰衣足食了.之前在学习法向量的计算时曾经说过,PCA算法最初其 ...
Open3D 点云曲率计算（Python版本）
文章目录一.简介二.代码实现三.实现效果参考资料一.简介由于在Open3D中也没发现相关函数,就自己动手丰衣足食了.之前在学习法向量的计算时曾经说过,PCA算法最初其实是一种评估曲率的方法 ...
点云曲率计算（MATLAB）
文章目录一.原理概述二.代码实现三.实现效果参考文献一.原理概述之前在学习法向量的计算时曾经说过,PCA算法最初其实是一种评估曲率的方法,当然表示一个点曲率的方法有很多种,但是一般他们都需 ...
Polygon Mesh Processing读书笔记——1三角网格Triangle Meshes
最近看论文深感基础知识的匮乏,所以补充一些图形几何方面的知识,首先是这本书的封面. 主要章节介绍本书讨论了基于多边形网格的几何处理管道的主要组件,如下图所示. 为了本书的指导目的,主题的描述顺序与图 ...
Siggraph三角网格变形之拉普拉斯变换
三角网格变形一直是CAGD相关领域的重点,刚上研究生的时候,感觉有点神奇.而且一上来导师就给我发了一篇基于格林坐标的自由变形的相关paper,让我看,外文文献,看了n多天,第一次看外文文献,啥也没看懂 ...
三角网格细分算法 —— Loop 算法
三角网格细分算法 -- Loop细分主要参考文章: https://zhuanlan.zhihu.com/p/144400261 https://blog.csdn.net/McQueen_LT/a ...
Delaunay 三角网格学习
本文是为<Mastering Opencv...>第七章准备的,他要使用主动外观模型AMM和POSIT算法做一个头部3D姿态估计.图像上的特征点由人工标定,但由于头部运动,比如张嘴,会导致 ...

三角网格的顶点曲率计算（平均曲率和高斯曲率）

参考：

曲率计算过程：

三角网格的顶点曲率计算（平均曲率和高斯曲率）相关推荐

最新文章

热门文章