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

  方法一，参考张学东的论文《空间曲线的曲率计算方法》：

  方法二，参考博文：空间曲线的弧长与曲率

二、空间曲线的曲率计算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('曲率半径')

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

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

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

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

参考：

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)