球面谐波函数 （Spherical harmonic function）分析实际颗粒形状公式推导及数值实现 Part 4 SH展开系数的重构

        使用球面谐波函数重构出颗粒的表面轮廓，并使用基于表面网格的形状特征分析方法，能够得到颗粒的形状特征，是球面谐波函数分析的正问题，当考虑到在离散元中使用不规则颗粒时，就需要研究如何借助SH展开的方法得到形状参数符合一定规律的颗粒这样的反问题。Wei等人的文献提供了一种简单的重构方法，并给出了相应的程序，

文章的链接如下：https://www.sciencedirect.com/science/article/abs/pii/S0032591018301189

代码的Github链接如下：https://github.com/budizhao/SHPSG

一、原理方法

二、球面描述符计算及程序

从系数矩阵中读取，并组合成

###        球谐系数读取与组合
### 作者：Chunqiwang 时间：2023-09-13
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.special import sph_harm
from mpl_toolkits.mplot3d import Axes3D
import openpyxl

## 一、读取系数，并将系数按照x,y,z三个方向分解
# 1. 打开文本文件并读取内容
with open('i_thParticle15.txt', 'r') as file:
    data_str = file.read()

# 2. 使用正则表达式提取所需的数据部分
import re

# 使用正则表达式查找匹配的内容
pattern = r'{([^{}]+)}'
matches = re.findall(pattern, data_str)

# 3. 将提取的数据转换为适当的数据结构
cx = []
cy = []
cz = []

for match in matches:
    elements = [float(x) for x in match.split(',')]
    cx.append(elements[0])
    cy.append(elements[1])
    cz.append(elements[2])

## 二、还原SH图像
'''
# 判断阶数
n = int(math.sqrt(np.size(cx)) - 1)
print(n)
# 创建一个网格点集合
theta = np.linspace(0, 2 * np.pi, 200)
phi   = np.linspace(0, np.pi, 100)
theta, phi = np.meshgrid(theta, phi)

xyz_2d = np.array([np.sin(phi) * np.sin(theta),
                   np.sin(phi) * np.cos(theta),
                   np.cos(phi)])
# 初始化还原的函数
reconstructed_functionx = np.complex64(np.zeros(theta.shape))
reconstructed_functiony = np.complex64(np.zeros(theta.shape))
reconstructed_functionz = np.complex64(np.zeros(theta.shape))

# x y z 分别计算
for i in range(n + 1):
    for j in range(-i,i + 1):
        N = i * (i + 1) + j
        yij = sph_harm(abs(j),i,theta,phi)
        reconstructed_functionx += cx[N] * yij
        reconstructed_functiony += cy[N] * yij
        reconstructed_functionz += cz[N] * yij
# 创建图形
print(reconstructed_functionz)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(reconstructed_functionx.real, reconstructed_functiony.real, reconstructed_functionz.real, rstride=1, cstride=1, antialiased=True, alpha=0.5)

plt.show()
'''

## 三、计算系数组合
dl = np.zeros(15)
for i in range(np.size(dl)):
    N = i + 1   # 阶数
    dlx = 0
    dly = 0
    dlz = 0
    for j in range(-N,N+1):
        index = N * (N + 1) + j
        dlx = dlx + np.abs(cx[index])   # 计算球面描述符的x分量
        dly = dly + np.abs(cy[index])
        dlz = dlz + np.abs(cz[index])
    dl[i] = dlx + dly + dlz
D2_8 = np.sum(dl[1:8]) / dl[0]
D9_15 = np.sum(dl[8:15]) / dl[0]

# 写入表格文件
# 创建一个新的工作簿和工作表
workbook = openpyxl.Workbook()
worksheet = workbook.active

# 将数据添加到表格的最后一行
worksheet.append(dl.tolist())

# 替换成你要保存的文件路径
excel_file_path = 'output.xlsx'   # 修改文件名
workbook.save(excel_file_path)

分析一个颗粒的和间的关系：

和Wei等人的文中类似，2-8和9-15的dl分别用一个函数进行拟合，得到系数α和β，并用于后续的计算。

三、SH系数重构

通过给定EI、FI、以及和通过上面拟合的式子可以得到每个dl，再通过系数构造方法，得到颗粒的轮廓。

1. 不同的EI（FI）参数和参数对颗粒形状的影响

一阶椭球的Ei和Fi仅影响颗粒的纵横比，不会对颗粒的圆度、粗糙度造成影响，D2_8会对颗粒的表面细微形状产生影响，且数值越大，表面越复杂，圆度越小。

2. 不同的对颗粒形状的影响

D9_15会对颗粒的表面粗糙度产生影响，不会对颗粒的纵横比产生影响，且数值越大，颗粒的粗糙程度越大。

3. 统计分析颗粒的形状参数与上述EI、FI、以及和的关系