球面谐波函数（Spherical harmonic function）分析实际颗粒形状公式推导及数值实现 Part 3 颗粒形状分析

老王博的科研笔记

已于 2023-09-12 22:59:42 修改

阅读量222

点赞数

文章标签： python 经验分享

于 2023-09-12 22:56:59 首次发布

本文链接：https://blog.csdn.net/weixin_41797057/article/details/132763324

版权

参考文献：3D quantitative shape analysis on form, roundness, and compactness
with μCT https://www.sciencedirect.com/science/article/pii/S0032591015302345

一、形状分析的种类

1. 宏观形状（Aspect ratio，AR 纵横比）

将颗粒的体素转置， $p_1,p_2,p_3$ 分别对应 $x,y,z$ 方向， $EI=p_2/p_1$ ， $FI=p_3/p_2$ ， $AR$ 为两者的平均值。

转置代码：

% 转置成x,y,z依次减小的格式
originalSize = size(i_th_Particle_Voxel);

% 将维度信息排序，以便它们依次减小
sortedSize = sort(originalSize, 'descend');
order = [1,1,1];
for i = 1:3
   for j = 1:3
        if sortedSize(i) == originalSize(j)
            order(i) = j;
        end
   end
end

i_th_Particle_Voxel = permute(i_th_Particle_Voxel, order);

2. 圆度

更为精细特征的形状描述符。使用基于曲面三角网格的局部曲率进行计算。具体方法是首先计算颗粒的内切球半径及曲率；然后计算每一个网格点的平均曲率，并用内切球曲率除以平均曲率得到局部圆度；随后循环每个三角面，将其三个顶点的局部圆度平均，并判断其是否>0&<1，只有满足这个条件，才能用与后续计算圆度；最后对满足条件的网格上的局部圆度对网格面积进行平均。最终得到圆度， $R_i$ ：

$\begin{equation} R_i=\frac{\sum\left(A_n \frac{k_s}{k_i}\right)}{\sum\left(A_n\right)} \end{equation}$

其中， $k_s$ 表示内切球曲率； $k_i$ 表示每个面的局部圆度； $A_n$ 表示每个三角面的面积。

3. 球度和密实度（compactness）

球度：V表示颗粒体积；SA表示颗粒轮廓表面积。

紧实度： $V_{CH}$ 表示包围颗粒最小凸包（convex hull）的体积

二、程序及结果展示

1. 程序框架

使用python的vtk库，对三角网格数据进行读入和计算其几何特征；安装vtk库：

pip install vtk

2. 程序代码（未整理版）

import vtk
import numpy as np
import math
from scipy.spatial import ConvexHull

# 创建一个vtk读取器
reader = vtk.vtkPolyDataReader()

# 设置vtk路径
vtk_file_path = "i_thParticle15.vtk"

# 将vtk文件加载到读取器中
reader.SetFileName(vtk_file_path)

# 执行读取操作
reader.Update()

# 读取vtk数据
vtk_data = reader.GetOutput()

# 创建一个vtk体积计算器
volume_calculator = vtk.vtkMassProperties()
volume_calculator.SetInputData(vtk_data)

# 执行体积计算
volume_calculator.Update()

# 获取体积
volume = volume_calculator.GetVolume()

# 打印体积
print('体积：')
print(volume)

#### 计算平均曲率，并进行可视化展示

centerOfMassFilter = vtk.vtkCenterOfMass()
centerOfMassFilter.SetInputData(vtk_data)
centerOfMassFilter.Update()

# 获取质心坐标
center = centerOfMassFilter.GetCenter()

# 获取曲面上的点
points = vtk_data.GetPoints()
num_points = points.GetNumberOfPoints()

# 初始化最小距离
min_distance = np.inf

# 计算每个点到质心的距离并找到最小距离
for i in range(num_points):
    point = np.array(points.GetPoint(i))
    distance = np.linalg.norm(point - center)
    if distance < min_distance:
        min_distance = distance

# 最小距离即为内切球半径
inscribed_sphere_radius = min_distance
print(f'内切球半径: {inscribed_sphere_radius}')

'''
# 创建一个vtk曲率计算器
curvature_calculator = vtk.vtkCurvatures()
curvature_calculator.SetInputData(vtk_data)
curvature_calculator.SetCurvatureTypeToGaussian()  # 设置曲率类型为高斯曲率
curvature_calculator.Update()

# 获取曲率数据
curvature_data = curvature_calculator.GetOutput().GetPointData().GetScalars()

# 创建一个空的NumPy数组来存储曲率数据
num_points = vtk_data.GetNumberOfPoints()
curvature_array = np.zeros(num_points)

# 将VTK数据中的曲率数据复制到NumPy数组中
for i in range(num_points):
    curvature_array[i] = curvature_data.GetValue(i)

# 指定曲率阈值，只保留曲率大于该阈值的点
curvature_threshold = 1 / inscribed_sphere_radius  # 替换为您的曲率阈值

# 筛选出满足曲率阈值条件的点
filtered_points = [vtk_data.GetPoint(i) for i in range(num_points) if curvature_array[i] > curvature_threshold]

# 计算满足条件点的平均曲率半径
average_curvature_radius = np.mean([1.0 / curvature_array[i] for i in range(num_points) if curvature_array[i] > curvature_threshold])

# 打印或输出满足条件点的平均曲率半径
print(f'满足曲率阈值条件的点的平均曲率半径: {average_curvature_radius}')
'''

##### 计算AR纵横比
# 获取点数据和面数据
points = vtk_data.GetPoints()
cells = vtk_data.GetPolys()
# 计算纵横比
min_point = np.min(points.GetData(), axis=0)
max_point = np.max(points.GetData(), axis=0)
scale = max_point - min_point
EI = scale[1] / scale[0]
FI = scale[2] / scale[1]
AR = (EI + FI) / 2.0
print('AR:')
print(AR)


##### 计算球度
# 创建一个vtk三角面片过滤器，用于分解曲面为三角面片
triangulate = vtk.vtkTriangleFilter()
triangulate.SetInputData(vtk_data)
triangulate.Update()

# 获取三角面片数据
triangles = triangulate.GetOutput()

# 计算曲面的表面积
surface_area = 0.0

for i in range(triangles.GetNumberOfCells()):
    cell = triangles.GetCell(i)
    area = cell.ComputeArea()
    surface_area += area

# 打印或输出表面积
print(f'曲面的表面积: {surface_area}')

# 计算球度
S = ((36 * math.pi * volume ** 2) ** (1/3)) / surface_area

print('球度：')
print(S)

### 密实度

#计算凸包体积
# 获取顶点坐标
points = vtk_data.GetPoints()
num_points = points.GetNumberOfPoints()

# 将顶点坐标转换为NumPy数组
point_array = np.zeros((num_points, 3))  # 创建一个NumPy数组来存储顶点坐标

for i in range(num_points):
    point = points.GetPoint(i)
    point_array[i] = point

# 使用scipy的ConvexHull计算凸包
hull = ConvexHull(point_array)

# 计算凸包体积
convex_hull_volume = hull.volume

# 打印或输出凸包体积
print(f'凸包的体积: {convex_hull_volume}')

# 计算紧实度
Cx = volume / convex_hull_volume
print('紧实度：')
print(Cx)



# 内切球曲率
ks = 1.0 / inscribed_sphere_radius

# 创建一个vtk曲率计算器
curvature_calculator = vtk.vtkCurvatures()
curvature_calculator.SetInputData(vtk_data)
curvature_calculator.SetCurvatureTypeToMean()  # 设置曲率类型为高斯曲率
curvature_calculator.Update()

# 获取曲率数据
curvature_data = curvature_calculator.GetOutput().GetPointData().GetScalars()

# 获取三角形网格的单元
cells = vtk.vtkCellArray()
cells = vtk_data.GetPolys()  # 对于三角形网格，使用GetPolys函数获取单元

# 创建一个数组来存储每个三角形中三个顶点的曲率平均值
triangle_vertex_curvature_averages = []

# 创建一个数组来存储每个三角形的面积
triangle_areas = []

# 遍历三角形网格的单元
cells.InitTraversal()
while True:
    cell = vtk.vtkIdList()
    if cells.GetNextCell(cell) == 0:
        break

    # 获取三角形的顶点坐标和曲率值
    triangle_curvatures1 = [curvature_data.GetValue(cell.GetId(i)) for i in range(cell.GetNumberOfIds())]
    triangle_curvatures = []
    for i in triangle_curvatures1:
        if i != 0:
            result = ks / i
            triangle_curvatures.append(result)
        else:
            triangle_curvatures.append(None)
    # 计算当前三角形中三个顶点的曲率平均值
    vertex_curvature_average = np.mean(triangle_curvatures)
    if (vertex_curvature_average >0) and (vertex_curvature_average < 1):
        # 将每个三角形中三个顶点的曲率平均值添加到结果数组中
        triangle_vertex_curvature_averages.append(np.abs(vertex_curvature_average))

        # 计算当前三角形的面积
        v1 = np.array(vtk_data.GetPoint(cell.GetId(0)))
        v2 = np.array(vtk_data.GetPoint(cell.GetId(1)))
        v3 = np.array(vtk_data.GetPoint(cell.GetId(2)))
        area = 0.5 * np.linalg.norm(np.cross(v2 - v1, v3 - v1))
        triangle_areas.append(area)

# 转换为NumPy数组
triangle_vertex_curvature_averages = np.array(triangle_vertex_curvature_averages)
# 计算总体加权平均曲率值
total_weighted_average_curvature = np.sum(triangle_vertex_curvature_averages * triangle_areas) / np.sum(triangle_areas)

# 打印总体加权平均曲率值
print(f'总体加权平均曲率值: {total_weighted_average_curvature}')

结果：

其中曲率的计算是本代码、复现的文章的重点，为此，第三小部分进行了可视化：

3. 可视化结果

# 可视化曲率

import vtk
import numpy as np
import math
from scipy.spatial import ConvexHull

# 创建一个vtk读取器
reader = vtk.vtkPolyDataReader()

# 设置vtk路径
vtk_file_path = "i_thParticle15.vtk"

# 将vtk文件加载到读取器中
reader.SetFileName(vtk_file_path)

# 执行读取操作
reader.Update()

# 读取vtk数据
vtk_data = reader.GetOutput()

# 创建一个vtk曲率计算器
curvature_calculator = vtk.vtkCurvatures()
curvature_calculator.SetInputData(vtk_data)
curvature_calculator.SetCurvatureTypeToMinimum()  # 设置曲率类型为高斯曲率 SetCurvatureTypeToMean 平均
# SetCurvatureTypeToGaussian 高斯  SetCurvatureTypeToMaximum 最大  SetCurvatureTypeToMinimum 最小
curvature_calculator.Update()

# 获取曲率数据
curvature_data = curvature_calculator.GetOutput().GetPointData().GetScalars()


# 将曲率数据添加到vtk数据中
vtk_data.GetPointData().SetScalars(curvature_data)

# 创建一个vtk写入器
writer = vtk.vtkPolyDataWriter()

# 设置输出vtk文件路径
output_file_path = "output_file_with_curvature.vtk"  # 替换为实际的输出文件路径
writer.SetFileName(output_file_path)

# 设置写入器的输入数据为带有曲率数据的vtk数据
writer.SetInputData(vtk_data)

# 执行写入操作
writer.Write()

每个点的平均曲率:

每个点的高斯曲率：

最大曲率：

最小曲率：