用python写一个元胞自动机四邻居的晶粒生长函数

好的，以下是用Python写一个元胞自动机四邻居的晶粒生长函数的代码：

import numpy as np
import matplotlib.pyplot as plt

# 初始化背景矩阵
# 0表示未生长的晶粒，1表示已生长的晶粒
# 初始状态设置为1 x 1的矩阵，即只有一个晶粒
background = np.array([[0]])

# 超参数，影响晶粒生长的速度和形态
prob = 0.3  # 生长概率
limit = 4   # 晶粒核心像素数量极限

# 生长函数
def growth(background):
    # 扩展背景，增加边缘
    bg_h, bg_w = background.shape
    background_expanded = np.zeros((bg_h+2, bg_w+2))
    background_expanded[1:bg_h+1, 1:bg_w+1] = background

    # 初始化新状态
    new_background = np.copy(background)

    # 遍历每一个像素
    for i in range(bg_h):
        for j in range(bg_w):
            pixel = background_expanded[i+1, j+1]
            neighbors = background_expanded[i:i+3, j:j+3]

            if pixel == 0:
                # 如果当前像素标记为未生长状态
                # 统计四邻域中生长状态（值为1）像素的数量
                growed_neighbors = np.count_nonzero(neighbors) - 1  # 减去中心像素

                # 如果四邻域中有生长状态的像素，计算扩散概率
                if growed_neighbors > 0 and np.random.rand() < prob:
                    # 如果扩散概率通过，更新为生长状态
                    # 晶粒核心像素数量不能超过限制
                    if growed_neighbors < limit:
                        new_background[i, j] = 1

    return new_background

# 迭代生长，生成多幅背景矩阵，模拟不同时间的晶粒生长状态
bg_list = [background]
for i in range(20):
    bg_list.append(growth(bg_list[-1]))

# 可视化结果
fig, axs = plt.subplots(4, 5, figsize=(8, 8))
for i in range(4):
    for j in range(5):
        axs[i, j].imshow(bg_list[i*5+j], cmap='gray')
        axs[i, j].axis('off')
        axs[i, j].set_title(f'Time {i*5+j}')
plt.tight_layout()
plt.show()

注意这只是一个简单的生长函数，主要是为了演示元胞自动机的原理。实际应用中，生长函数的实现可能更加复杂，可能需要考虑更多因素的影响，如温度、压力等。