参考文献:3D quantitative shape analysis on form, roundness, and compactness
with μCT https://www.sciencedirect.com/science/article/pii/S0032591015302345
一、形状分析的种类
1. 宏观形状(Aspect ratio,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,只有满足这个条件,才能用与后续计算圆度;最后对满足条件的网格上的局部圆度对网格面积进行平均。最终得到圆度,:
其中,表示内切球曲率;
表示每个面的局部圆度;
表示每个三角面的面积。
3. 球度和密实度(compactness)
球度:V表示颗粒体积;SA表示颗粒轮廓表面积。
紧实度:表示包围颗粒最小凸包(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()
每个点的平均曲率:
每个点的高斯曲率:
最大曲率:
最小曲率: