今回は雪の結晶を作るプログラムです。

解説動画はこちら


ありのままの雪の結晶を作ってみたくて
雪の結晶の生成方法を探していたら
こんなサイトがありました。

参考:Pythonでフラクタル氷晶の
成長をシミュレートする方法
フーリエ変換と密度関数

http://articlesbyaphysicist.com/snowflake.html

要約すると

空気中に遊離していた水分子が氷の結晶に着地して付着すると
氷の結晶が成長する可能性があります。

結晶の表面にギザギザの山と深い谷があると想像すると
さまよう水分子は谷よりもはるかに山にぶつかる
可能性が高いため山はより速く成長します。

これにより、結晶が完全な球体として成長する可能性は低くなり
それらは周囲に棘やスパイクができてしまいます。

これがどのように機能するかをシミュレートすると
雪の結晶の氷の結晶に似たフラクタルが得られます。

だそうです。
何言っているのか分からないですねwwwww

ということでコードを拝借しまして
少し直させていただきました。

import numpy as np
import math
from PIL import Image

def red(x):
    return np.uint(math.exp(-x * 4.0) * 255)

def green(x):
    return np.uint(math.exp(-x) * 255)

def blue(x):
    return np.uint(math.exp(-x / 4.0) * 255)

def getN(n):
    rs = np.zeros([n, n])
    rs[0, 0],rs[0, 1],rs[1, 0],rs[-1, 0], rs[0, -1],rs[1, -1],rs[-1, 1] = 1,1,1,1,1,1,1
    return rs

def getG(n):
    G = np.zeros([n, n])
    for i in range(n):
        for j in range(n):
            i2 , j2 = i , j
            if i2 > n/2:
                i2 -= n
            if j2 > n/2:
                j2 -= n
            dx , dy = i2 + 0.5 * j2 , j2 * math.sqrt(3)/2
            r = math.sqrt(dx*dx + dy*dy)
            G[i, j] = math.exp(-r*r/8)
    return G / np.sum(np.sum(G))

def save_img(A , frame):
    B , n = A * 0.2 , len(A)
    B[B < 0.01] = 0
    B = np.array([B[i, range((n // 2 - i // 2), (n - i // 2), 1)] for i in range(n // 4, (3 * n) // 4, 1) if (i % 5) != 0])
    col = [[[red(x), green(x), blue(x), 255] for x in y] for y in B]
    im = Image.fromarray(np.uint8(col))
    im.save("frames/im{0:06}.png".format(frame))
    return frame + 1

def main(n=500 , rg=1000):
    A , frame = np.zeros([n, n]) , 0
    A[n//2, n//2] = 1
    G , nums = getG(n) , getN(n)
    rho = np.ones([n, n])
    sf , gf = np.fft.fft2(nums) , np.fft.fft2(G)
    for i in range(rg):
        B = A.copy()
        B[B < 1] = 0
        B[B > 1] = 1
        A2 = np.fft.ifft2(np.multiply(sf, np.fft.fft2(B))).real
        A2[A2 > 0.9] = 1
        A2[A2 < 0.9] = 0
        A += A2 * (rho * 0.5)
        A[A > 5] = 5
        rho[A > 0.5] = 0
        rho = np.fft.ifft2(np.multiply(gf, np.fft.fft2(rho))).real
        rho = rho / np.mean(np.mean(rho)) + np.random.randn(n, n) * 0.02
        if (i%10)==0:
            frame = save_img(A,frame)
    frame = save_img(A,frame)

先に「frames」というフォルダを
プログラム配下に用意してください。

実行は画像サイズと
繰り返しの回数を指定するだけです。

# 要 frames フォルダ
size = 800
nums = 2000
main(size,nums)

結果は「frames」フォルダ内に
画像ファイルが出力されます。

画像を見てみましょう。
ウィジェットで好きな項番を見られます。
from ipywidgets import interact,IntSlider
import matplotlib.pyplot as plt
%matplotlib inline

frame_num = IntSlider(min=0, max=200, step=1, value=0)
@interact(frame_num=frame_num)
def plot(frame_num):
    fig = plt.figure(figsize=(10, 10))
    img = plt.imread("frames/im{0:06}.png".format(frame_num))
    plt.imshow(img)
    plt.axis('off')
    plt.show()
スクリーンショット 2021-11-27 17.18.21



最後にこれをGIF画像にしてみましょう。
間引いてGIFにします。
from PIL import Image
import os
import glob

out_file_path='snowflake.gif'
imgs = []
for file_path in sorted(glob.glob(os.path.join(*['frames', '*.png']))) :
    img = Image.open(file_path)
    imgs.append(img) 
imgs[0].save(out_file_path,save_all=True, append_images=imgs[1::10], optimize=False, duration=300, loop=0)

出来上がったGIFを見てみましょう。

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
from IPython import display
from pathlib import Path

file_path = Path("snowflake.gif")
with open(file_path,'rb') as f:
    display.Image(data=f.read(), format='png')
snowflake