Hello World / plɹoM ollǝH

Programmers Live in Vain

Turing Patterns

Gray-Scottモデルの反応拡散式で模様を自動生成してみました。

f:id:dungeonneko:20150624223537p:plain

C++

pythonでは処理速度が厳しいため反応拡散の計算はC++で書きました。

#define BOOST_PYTHON_STATIC_LIB
#include <boost/python.hpp>
#include <vector>
#include <algorithm>

boost::python::list Create(int in_w, int in_h, float in_du, float in_dv, float in_f, float in_k, int in_t)
{
    float* bufferU_0 = new float[in_w * in_h];
    float* bufferU_1 = new float[in_w * in_h];
    float* bufferV_0 = new float[in_w * in_h];
    float* bufferV_1 = new float[in_w * in_h];

    float* bufferU[2];
    float* bufferV[2];

    bufferU[0] = bufferU_0;
    bufferU[1] = bufferU_1;
    bufferV[0] = bufferV_0;
    bufferV[1] = bufferV_1;

    int current_idx = 0;

    // Initialize
    for (int y = 0; y < in_h; ++y) {
        for (int x = 0; x < in_w; ++x) {
            int c = y * in_w + x;
            bufferU[current_idx][c] = 1.0f;
            bufferV[current_idx][c] = (rand() % 100) < 3 ? 1.0f : 0.0f;
        }
    }

    // Compute
    for (int t = 0; t < in_t; ++t) {
        int next_idx = (current_idx + 1) % 2;
        for (int y = 0; y < in_h; ++y) {
            for (int x = 0; x < in_w; ++x) {
                int u = std::max(y - 1, 0) * in_w + x;
                int d = std::min(y + 1, in_h - 1) * in_w + x;
                int l = y * in_w + std::max(x - 1, 0);
                int r = y * in_w + std::min(x + 1, in_w - 1);
                int c = y * in_w + x;

                float uu = bufferU[current_idx][u];
                float ud = bufferU[current_idx][d];
                float ul = bufferU[current_idx][l];
                float ur = bufferU[current_idx][r];
                float uc = bufferU[current_idx][c];
                float vu = bufferV[current_idx][u];
                float vd = bufferV[current_idx][d];
                float vl = bufferV[current_idx][l];
                float vr = bufferV[current_idx][r];
                float vc = bufferV[current_idx][c];

                float diff_u = in_du * (uu + ud + ur + ul - 4.0f * uc);
                float diff_v = in_dv * (vu + vd + vr + vl - 4.0f * vc);
                float reac_x = uc * vc * vc;
                float reac_u = in_f * (1.0f - uc);
                float reac_v = (in_f + in_k) * vc;
                uc += diff_u - reac_x + reac_u;
                vc += diff_v + reac_x - reac_v;

                bufferU[next_idx][c] = std::min(1.0f, std::max(0.0f, uc));
                bufferV[next_idx][c] = std::min(1.0f, std::max(0.0f, vc));
            }
        }
        current_idx = next_idx;
    }

    // Finalize
    boost::python::list ret;
    for (int y = 0; y < in_h; ++y) {
        for (int x = 0; x < in_w; ++x) {
            ret.append(bufferV[current_idx][y * in_w + x]);
        }
    }

    delete[] bufferU_0;
    delete[] bufferU_1;
    delete[] bufferV_0;
    delete[] bufferV_1;

    return ret;
}

BOOST_PYTHON_MODULE(gs)
{
    boost::python::def("Create", &Create);
}

Python

PILを使用してC++の計算結果を画像データに変換しています。

# conding: utf-8
from PIL import Image
import gs

w = 128
h = 128
Du = 0.1
Dv = 0.05
f = 0.04
k = 0.06
t = 1000

img = Image.new('L', (w, h))
p = img.load()
srcImg = gs.Create(w, h, Du, Dv, f, k, t)
for y in range(h):
    for x in range(w):
        p[ x, y ] = int(srcImg[ y * w + x ] * 255)
img.show()

パラメータ設定が難しい。