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


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

$$ 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
```

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

![3_0.png](https://qiita-image-store.s3.amazonaws.com/0/147715/16753c65-5b2c-3864-3910-f9678d115639.png)
![3_1.png](https://qiita-image-store.s3.amazonaws.com/0/147715/3b59a6c6-261e-1df2-324f-4ef6deca209f.png)
![3_2.png](https://qiita-image-store.s3.amazonaws.com/0/147715/585c6cf4-3cae-85be-1fd6-9b04ffd6acb1.png)
![3_3.png](https://qiita-image-store.s3.amazonaws.com/0/147715/36e0344e-55cb-0127-f006-094b40ff0e47.png)
![3_4.png](https://qiita-image-store.s3.amazonaws.com/0/147715/a00549ee-601e-279c-b1c1-5ec682f3b660.png)
![3_5.png](https://qiita-image-store.s3.amazonaws.com/0/147715/b412fafd-effa-8799-098e-f8a1c5649096.png)
![3_6.png](https://qiita-image-store.s3.amazonaws.com/0/147715/cb264f8c-2410-5c6c-bffc-2c9c37a46b44.png)
![3_7.png](https://qiita-image-store.s3.amazonaws.com/0/147715/834649f0-bf3a-a72a-1752-2630e8e8a8dc.png)
![3_8.png](https://qiita-image-store.s3.amazonaws.com/0/147715/89bd6648-2c3e-329b-3ec0-ebd0ed8293a2.png)
![3_9.png](https://qiita-image-store.s3.amazonaws.com/0/147715/03c37c15-d013-2cd3-3a4d-7834867d2d37.png)
![3_10.png](https://qiita-image-store.s3.amazonaws.com/0/147715/4bed2b83-e7f1-192d-46e4-791c705fe74a.png)