如何从3D温度numpy数组中提取T=900等值面索引以获取对应梯度值

你可以继续使用已经引入的skimage.measure模块，调用其内置的3D等值面提取接口marching_cubes实现需求，无需额外引入绘图类依赖，完整实现逻辑如下：

import numpy as np
from skimage import measure

# 你原有生成网格、温度场、梯度的代码保持不变
res = 0.1 # 按你实际的分辨率赋值即可
def regular(x,y,z,q=100,k=175,a=7.1e-5):
    R = np.sqrt(x**2+y**2+z**2)
    return 100 / (2*np.pi*k) * (1/R) * np.exp(-0.5/a*(R+x))
x = np.arange(-1.5,0.5+res/2,res)*1e-3
y = np.arange(-1.0,1.0+res/2,res)*1e-3
z = np.arange(0.0,0.5+res/2,res)*1e-3
Y,X,Z = np.meshgrid(y,x,z)
T = regular(X,Y,Z)
dTdx, dTdy, dTdz = np.gradient(T) # 修正原伪代码中的dYdz笔误

# 提取T=900的等值面
verts, faces, normals, values = measure.marching_cubes(T, level=900)
# verts为形状(N,3)的数组，每一行对应等值面上一个点的坐标，默认顺序为(z,y,x)

# 若需要和你原有2D方案逻辑一致，取最近的整数网格索引，直接四舍五入转整型即可
verts_int = np.round(verts).astype(int)
zind, yind, xind = verts_int[:,0], verts_int[:,1], verts_int[:,2]

# 直接索引得到等值面对应位置的梯度值
x_gradients_at_isosurface = dTdx[xind, yind, zind]
y_gradients_at_isosurface = dTdy[xind, yind, zind]
z_gradients_at_isosurface = dTdz[xind, yind, zind]

注意事项：

• measure.marching_cubes返回的顶点坐标默认顺序为(z, y, x)，需和你生成温度场T时的维度顺序匹配，若你实际T的维度顺序不同，调整索引取值顺序即可
• 若需要更高精度的等值面梯度值，不需要四舍五入取整，可以直接用verts返回的浮点坐标对三个梯度数组做三线性插值
• 该接口会自动返回所有不连通的T=900等值面的坐标，无需额外处理多连通区域

内容的提问来源于stack exchange，提问作者F. Heath