水面に垂らしたインクの拡散シミュレーション

以下の拡散方程式を用いて作りました。

$$ U_{i,j}^{t+1} - U_{i,j}^{t} = D \times ( U_{i-1,j}^t - 2U_{i,j}^t + U^t_{i+1,j} ) + D \times ( U^t_{i,j-1} - 2U^t_{i,j} + U^t_{i,j+1})$$

パストレーシングに手を出そうとしているので練習を兼ねて、ppm形式の画像ファイルを出力しています。

InkFlow.java

package flow;

import java.io.FileWriter;
import java.io.IOException;

public class InkFlow{
    // 水面に落としたインクの拡散をシミュレーションする
    static double INK = 255.; // RGBの最大値
    static int mx; // 横のメッシュ数(画像の画素数？になる)
    static int my; // 縦
    static int fc=0; // ファイル名の変数

    public static void main(String[] args) throws IOException {
        /*******
         * 設定(ここを変える)
         *******/
        // 水面の設定
        int x = 10; // 横方向の長さ[cm]
        int y = 10; // 縦方向の長さ[cm]

        // インクを垂らす場所（左上、上、右上、左、中央、右、左下、下、右下の９箇所）
        boolean CPlace[][] = { // Cyan
                {false,false,false},
                {true,false,false},
                {false,false,false}
        };
        boolean MPlace[][] = { // Magenta
                {false,true,false},
                {false,false,false},
                {false,false,false}
        };
        boolean YPlace[][] = { // Yellow
                {false,false,false},
                {false,false,false},
                {false,true,false}
        };
        // 垂らすインクの大きさ
        int InkSize = 13;

        // 時間の設定
        int t=0; // 初期時刻0[msec]
        int T=10000;     // 終了時刻[msec]
        int ut = 1000;      // ファイル出力の単位時間[msec]
        int inkt = 5000; // インクを垂らしている時間[msec]

        // 拡散用の変数
        double ul = 0.1; // 単位距離[cm]
        mx = (int) (x / ul);
        my = (int) (y / ul);
        double D = 0.25; // 拡散係数(<0.25)

        LogSettings(x,y,InkSize,T,ut, inkt ,ul,mx,my,D); // 設定のログ

        /*******
         * 初期化
         *******/
        // 色の初期化（白）
        double R[][] = InitialInk();
        double G[][] = InitialInk();
        double B[][] = InitialInk();

        // インクを垂らす
        DipInk(CPlace, MPlace, YPlace, InkSize, R, G, B);
        Output(R, G, B);

        /*******
         * メインループ
         *******/
        while(t++ <= T) { // 終了時刻Tになるまでループ
            // 拡散
            R = CalcDiffusion(R, D);
            G = CalcDiffusion(G, D);
            B = CalcDiffusion(B, D);
            // インクを垂らす
            if(t < inkt) { // inkt msec垂らす
                DipInk(CPlace, MPlace, YPlace, InkSize, R, G, B);
            }
            // 出力
            if(t % ut == 0){
                Output(R, G, B);
            }
        }
        System.out.println("done");
    }

    /*******
     * 各種メソッド
     *******/
    static double[][] InitialInk(){ // インクの初期化
        double a[][] = new double[mx][my];
        for(int i = 0; i < mx; i++) {
            for(int j = 0; j < my; j++) {
                a[i][j] = INK;
            }
        }
        return a;
    }

    static double[][] CalcDiffusion(double[][] a, double D){ // 拡散の計算
        double b[][] = new double[mx][my];
        // 内側の計算
        for(int i=1;i<(mx-1);i++){
            for(int j=1;j<(my-1);j++){ // 拡散方程式
                b[i][j] = a[i][j] + D * (a[i+1][j] + a[i-1][j] - 2*a[i][j]) + D * (a[i][j+1] + a[i][j-1] - 2*a[i][j]);
            }
        }
        // 外側の計算
        for(int i=0;i<mx;i++){
            b[i][my-1] = b[i][my-2];
            b[i][0] = b[i][1];
        }
        for(int j=0;j<my;j++){
            b[mx-1][j] = b[mx-2][j];
            b[0][j] = b[1][j];
        }
        return b;
    }

    static void DipInk(boolean[][] CPlace, boolean[][] MPlace, boolean[][] YPlace, int InkSize, double[][] R, double[][] G, double[][] B) {
        for(int i = 0; i < CPlace.length; i++) {
            for(int j = 0; j < CPlace[0].length; j++) {
                if(CPlace[i][j] || MPlace[i][j] || YPlace[i][j]) {
                    if(CPlace[i][j]) {
                        DipWithInkSize(R, InkSize, i, j, 0);
                    } else {
                        DipWithInkSize(R, InkSize, i, j, 212);
                    }
                    if(MPlace[i][j]) {
                        DipWithInkSize(G, InkSize, i, j, 0);
                    } else {
                        DipWithInkSize(G, InkSize, i, j, 212);
                    }
                    if(YPlace[i][j]) {
                        DipWithInkSize(B, InkSize, i, j, 0);
                    } else {
                        DipWithInkSize(B, InkSize, i, j, 212);
                    }
                }
            }
        }
    }
    static void DipWithInkSize(double[][] C, int InkSize, int a, int b, double val) {
        // インクサイズに応じて、インクのマスを渦巻き状に増やす
        C[(int) (mx*(2*a+1)/6)][(int) (my*(2*b+1)/6)] = val; // 1つ目のインク
        int itr = 1;
        int i = 1;
        while(true) {
            for(int x = itr, y = 0; x > 0; x--,y++) {
                C[(int) (mx*(2*a+1)/6)+x][(int) (my*(2*b+1)/6)+y] = val;
                i++;
                if(i >= InkSize) {
                    return;
                }
            }
            for(int x = 0, y = itr; y > 0; x--, y--) {
                C[(int) (mx*(2*a+1)/6)+x][(int) (my*(2*b+1)/6)+y] = val;
                i++;
                if(i >= InkSize) {
                    return;
                }
            }
            for(int x = -itr, y = 0; x < 0; x++,y--) {
                C[(int) (mx*(2*a+1)/6)+x][(int) (my*(2*b+1)/6)+y] = val;
                i++;
                if(i >= InkSize) {
                    return;
                }
            }
            for(int x = 0, y = -itr; y < 0; x++, y++) {
                C[(int) (mx*(2*a+1)/6)+x][(int) (my*(2*b+1)/6)+y] = val;
                i++;
                if(i >= InkSize) {
                    return;
                }
            }
            itr++;
        }
    }

    static void Output(double[][] R, double[][] G, double[][] B) throws IOException { // 画像ファイル(ppm形式)の出力
        String fname = "3_" + fc++ + ".ppm";
        FileWriter fw = new FileWriter(fname);    // Difine file name
        // Header
        fw.write("P3 \r\n");
        fw.write("#The P3 means colors are in ASCII, then columns and rows, then 255 for max color, then RGB triplets \r\n");
        fw.write(mx + " " + my + " \r\n");
        fw.write("" + (int) INK + " \r\n");

        // Body
        for(int i = 0; i < mx; i++) {
            for(int j = 0; j < my; j++) {
                fw.write("" + (int) R[i][j] + " " + (int) G[i][j] + " " + (int) B[i][j]);
                fw.write(" \r\n");
            }
        }
        fw.close();
        System.out.println("output: " + fname);
    }

    static void LogSettings(int x, int y, int InkSize, int T, int ut, int inkt, double ul, int mx, int my, double D) {
        System.out.println("水面(横×縦)  : " + x + " cm × " + y + " cm");
        System.out.println("インクサイズ  : " + ul*ul*InkSize + " cm^2");
        System.out.println("垂らす時間    :" + inkt + " msec");
        System.out.println("終了時間     : " + T + " msec");
        System.out.println("出力単位時間  : " + ut + " msec");
        System.out.println("単位距離     : " + ul + " cm");
        System.out.println("横のメッシュ数: " + mx);
        System.out.println("縦のメッシュ数: " + my);
        System.out.println("拡散係数     : " + D);
        System.out.println("-*-*-*-");
        if(mx / 6 < InkSize) {
            System.out.println("InkSizeが大きすぎます");
        }
    }
}

consoleでは以下のように出力されます。

水面(横×縦)  : 10 cm × 10 cm
インクサイズ  : 0.13000000000000003 cm^2
垂らす時間    :5000 msec
終了時間     : 10000 msec
出力単位時間  : 1000 msec
単位距離     : 0.1 cm
横のメッシュ数: 100
縦のメッシュ数: 100
拡散係数     : 0.25
-*-*-*-
output: 3_0.ppm
output: 3_1.ppm
output: 3_2.ppm
output: 3_3.ppm
output: 3_4.ppm
output: 3_5.ppm
output: 3_6.ppm
output: 3_7.ppm
output: 3_8.ppm
output: 3_9.ppm
output: 3_10.ppm
done

出力される画像は以下のようなものに