曲率与法线

高斯曲率（Gaussian Curvature）

定义

高斯曲率K 是通过两个主曲率（主方向上的曲率）的乘积来定义的。设曲面S 上的某点p 处有两个正交的主曲率k_1 和k_2 ，则该点的高斯曲率定义为：

K = k_1 \cdot k_2

性质

• 正高斯曲率：如果K > 0 ，曲面在该点呈现椭圆形，如球面上的点。
• 零高斯曲率：如果K = 0 ，曲面在该点是平坦的，如圆柱面或平面上的点。
• 负高斯曲率：如果K < 0 ，曲面在该点呈现鞍形，如马鞍面上的点。

例子

• 球面：在球面上任意一点，两个主曲率是相等的，且均为正值，因此高斯曲率总是正的。
  对于半径为r 的球面，任意一点的高斯曲率为：
  K = \frac{1}{r^2}
• 平面：平面上的任意一点，两个主曲率均为零，因此高斯曲率也是零。
• 圆柱面：圆柱面上的一点，其主曲率之一为零，另一个为非零常数，因此高斯曲率为零。
• 鞍面：在鞍面上，两个主曲率分别为正和负，因此高斯曲率为负。

高斯曲率的重要性

高斯曲率是内在几何量，它只依赖于曲面的第一基本形式，与曲面的嵌入方式无关。高斯的标志性定理（Theorema Egregium）表明，高斯曲率是一个内在属性，这意味着它可以通过曲面的本身来测量，而不依赖于曲面在空间中的具体位置。

平均曲率（Mean Curvature）

定义

平均曲率H 是两个主曲率k_1 和k_2 的算术平均值。设曲面S 上的某点p 处有两个主曲率k_1 和k_2 ，则该点的平均曲率定义为：

H = \frac{k_1 + k_2}{2}

性质

• 平均曲率可以是正、负或零，这取决于两个主曲率的符号和大小。
• 当平均曲率为零时，曲面称为极小曲面（minimal surface）。

例子

• 球面：在球面上任意一点，两个主曲率是相等的，因此平均曲率为：
  H = \frac{k_1 + k_2}{2} = \frac{2 \cdot \frac{1}{r}}{2} = \frac{1}{r}
• 平面：平面上的任意一点，两个主曲率均为零，因此平均曲率也是零。
• 圆柱面：圆柱面上的一点，其主曲率之一为零，另一个为非零常数k ，因此平均曲率为：
  H = \frac{0 + k}{2} = \frac{k}{2}
• 鞍面：在鞍面上，两个主曲率分别为正和负，如果k_1 = -k_2 ，则平均曲率为零。

平均曲率的重要性

平均曲率在物理学和几何学中有广泛应用。例如，极小曲面（如肥皂膜）是表面张力最小的曲面，其平均曲率为零。平均曲率还用于研究曲面的稳定性和形变问题。

求高斯曲率和法线，并各自绘图

import numpy as np
import open3d as o3d
import matplotlib.pyplot as plt

# Generate a point cloud
x = np.linspace(-2, 2, 400)
y = np.linspace(-2, 2, 400)
x, y = np.meshgrid(x, y)
z = np.sin(np.sqrt(x ** 2 + y ** 2))

xyz = np.zeros((np.size(x), 3))
xyz[:, 0] = np.reshape(x, -1)
xyz[:, 1] = np.reshape(y, -1)
xyz[:, 2] = np.reshape(z, -1)

pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(xyz)

# Estimate normals
pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30))

# Compute curvature for each point using PCA
def compute_curvature(pcd):
    points = np.asarray(pcd.points)
    normals = np.asarray(pcd.normals)
    curvatures = np.zeros(len(points))
    
    kdtree = o3d.geometry.KDTreeFlann(pcd)
    
    for i in range(len(points)):
        # Get the neighbors of the point
        [k, idx, _] = kdtree.search_knn_vector_3d(points[i], 30)
        neighbors = points[idx[1:], :]  # Exclude the point itself
        
        # Compute covariance matrix
        cov = np.cov(neighbors - points[i], rowvar=False)
        
        # Compute eigenvalues
        eigenvalues, _ = np.linalg.eigh(cov)
        
        # Curvature is the smallest eigenvalue
        curvatures[i] = eigenvalues[0] / np.sum(eigenvalues)
    
    return curvatures

curvatures = compute_curvature(pcd)

# Print the normal and curvature of the first point
print("Normal of the first point: ", pcd.normals[0])
print("Curvature of the first point: ", curvatures[0])

# Visualize the point cloud with curvature-based coloring
curvature_colors = plt.get_cmap('plasma')(curvatures / np.max(curvatures))
curvature_colors = curvature_colors[:, :3]  # Get RGB values

pcd.colors = o3d.utility.Vector3dVector(curvature_colors)

# Create a visualization window and draw the normals
o3d.visualization.draw_geometries([pcd], window_name="Point Cloud with Curvature", width=800, height=600)

# Draw normals
normals = pcd.normals
points = pcd.points

line_set = o3d.geometry.LineSet()

# Starting points of the normals are the points of the point cloud
line_set.points = points

# Ending points of the normals are the points of the point cloud plus the normal vectors
lines = [[i, i + len(points)] for i in range(len(points))]
line_set.lines = o3d.utility.Vector2iVector(lines)

colors = [[0, 1, 0] for i in range(len(lines))]  # Color of the normals
line_set.colors = o3d.utility.Vector3dVector(colors)

# Coordinates of the end points of the normals
end_points = points + np.asarray(normals) * 0.1

# Merge points and end points
all_points = np.vstack((points, end_points))
line_set.points = o3d.utility.Vector3dVector(all_points)

# Visualize the point cloud with normals
o3d.visualization.draw_geometries([pcd, line_set], window_name="Point Cloud with Curvature and Normals", width=800, height=600)