Python 四随机数生长法生成三维多孔介质代码及可视化

本代码使用四随机数生长法生成三维多孔介质，并利用 matplotlib 库将其可视化。

代码:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection


def generate_porous_medium(size, porosity):
    # 初始化
    medium = np.zeros(size)
    num_voxels = np.prod(size)
    num_voids = int(num_voxels * porosity)
    void_indices = np.random.choice(num_voxels, size=num_voids, replace=False)
    medium.flat[void_indices] = 1

    # 生长
    for i in range(10):
        # 随机打乱处理顺序
        indices = np.random.permutation(num_voxels)
        for index in indices:
            # 如果当前位置是空气孔隙，则跳过
            if medium.flat[index] == 1:
                continue

            # 获取当前位置的坐标
            x, y, z = np.unravel_index(index, size)

            # 随机选择一个方向
            direction = np.random.randint(6)
            dx, dy, dz = get_directions(direction)

            # 获取相邻位置的值
            x2 = x + dx
            y2 = y + dy
            z2 = z + dz
            if x2 < 0 or x2 >= size[0] or y2 < 0 or y2 >= size[1] or z2 < 0 or z2 >= size[2]:
                continue
            neighbor_value = medium[x2, y2, z2]

            # 根据相邻位置的值和概率确定当前位置的值
            if neighbor_value == 0:
                prob = 0.2
            else:
                prob = 0.8
            if np.random.rand() < prob:
                medium[x, y, z] = 1

    return medium


def get_directions(direction):
    if direction == 0:
        return -1, 0, 0
    elif direction == 1:
        return 1, 0, 0
    elif direction == 2:
        return 0, -1, 0
    elif direction == 3:
        return 0, 1, 0
    elif direction == 4:
        return 0, 0, -1
    elif direction == 5:
        return 0, 0, 1


def plot_porous_medium(medium):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    # 获取表面
    surface = get_surface(medium)

    # 绘制多边形集合
    collection = Poly3DCollection(surface, alpha=0.2)
    collection.set_facecolor('gray')
    ax.add_collection3d(collection)

    # 设置坐标轴范围
    ax.set_xlim(0, medium.shape[0])
    ax.set_ylim(0, medium.shape[1])
    ax.set_zlim(0, medium.shape[2])

    plt.show()


def get_surface(medium):
    surface = []

    # 遍历每个方格
    for x in range(medium.shape[0]):
        for y in range(medium.shape[1]):
            for z in range(medium.shape[2]):
                if medium[x, y, z] == 0:
                    # 如果当前方格是实体，则检查相邻的方格是否是空气孔隙
                    if x == 0 or medium[x - 1, y, z] == 1:
                        surface.append([(x, y, z + 1), (x + 1, y, z + 1), (x + 1, y + 1, z + 1), (x, y + 1, z + 1)])
                    if x == medium.shape[0] - 1 or medium[x + 1, y, z] == 1:
                        surface.append([(x + 1, y, z + 1), (x + 1, y + 1, z + 1), (x + 1, y + 1, z), (x + 1, y, z)])
                    if y == 0 or medium[x, y - 1, z] == 1:
                        surface.append([(x, y, z), (x + 1, y, z), (x + 1, y, z + 1), (x, y, z + 1)])
                    if y == medium.shape[1] - 1 or medium[x, y + 1, z] == 1:
                        surface.append([(x, y + 1, z + 1), (x + 1, y + 1, z + 1), (x + 1, y + 1, z), (x, y + 1, z)])
                    if z == 0 or medium[x, y, z - 1] == 1:
                        surface.append([(x, y, z), (x, y + 1, z), (x + 1, y + 1, z), (x + 1, y, z)])
                    if z == medium.shape[2] - 1 or medium[x, y, z + 1] == 1:
                        surface.append([(x + 1, y, z + 1), (x + 1, y + 1, z + 1), (x, y + 1, z + 1), (x, y, z + 1)])

    return surface


# 示例用法
size = (50, 50, 50)
porosity = 0.5
medium = generate_porous_medium(size, porosity)
plot_porous_medium(medium)

示例用法:

运行示例代码会生成一个大小为(50, 50, 50)、孔隙率为0.5的多孔介质，然后在一个新窗口中显示它。你可以尝试不同的参数来生成不同的多孔介质。

代码解释:

初始化: 代码首先创建一个大小为 size 的零矩阵，代表初始的多孔介质。然后根据 porosity 计算出需要生成的孔隙数量，并随机选择 num_voids 个索引作为孔隙。
生长: 循环执行 10 次生长过程。每次生长过程中，会遍历所有方格，并根据一定的概率规则确定每个方格是否成为孔隙。
可视化: 使用 matplotlib 库将生成的介质可视化。代码使用 Poly3DCollection 类绘制表面，并将表面颜色设置为灰色。

注意:

代码中 prob 的值影响生成的孔隙率。prob 越大，孔隙率越高。
生长次数（10 次）和 prob 的值可以根据实际需求调整。
此代码仅为示例，实际应用中可能需要根据具体情况进行修改。

本代码使用简单的四随机数生长法，可以快速生成三维多孔介质。但此方法生成的介质可能比较粗糙，难以模拟真实的复杂结构。对于需要更精细的多孔介质模型，建议使用更高级的算法。

更多信息: