EclipseでImageJのPlugin作成 -２次元拡散方程式シミュレーション-
前回は、モンテカルロ法による分子拡散のシミュレーションを書いたので
EclipseでImageJのPlugin作成 -拡散シミュレーション(ランダムウォーク)- - 生物屋さんのためのゼロからのプログラミング、今回は、２次元拡散方程式のシミュレーションを行ってみた。
(但し、普通の２次元拡散方程式のシミュレーションではなく、顕微鏡で観察した場合を想定したシミュレーションになっている。)
(注：メモリを非常に使うため、"java.lang.OutOfMemoryError: Java heap space"のエラーが出やすいので、Time等は短くしないとシミュレーションできない)
コードは下記。
import java.awt.BorderLayout;
import java.awt.Color;
import java.awt.Dimension;
import java.awt.Font;
import java.awt.Toolkit;
import java.awt.event.ActionEvent;
import java.awt.event.ActionListener;
import javax.swing.JButton;
import javax.swing.JFrame;
import javax.swing.JLabel;
import javax.swing.JPanel;
import javax.swing.JTextField;
import ij.IJ;
import ij.ImagePlus;
import ij.plugin.PlugIn;
import ij.process.ImageProcessor;
public class Fick_Diffusion implements PlugIn, ActionListener {
ImagePlus imp0, imp;
ImageProcessor ip;
JFrame frame;
JTextField dcField; //拡散係数
JTextField cellField; //カメラの画素サイズ
JTextField objField; //対物レンスの倍率
JTextField rateField; //撮影間隔
JTextField timeField; //撮影時間
int cellSize;
int obj;
int rate;
int localNumber;
int number;
int size;
float dt = 0.0005f; //0.5msec
float dc;
float pixelSize;
float speed;
float rawData [][][];
float data [][][];
////////////////////////////////////////////////
//キャンバスの作成
public void makeCanvas () {
//リセット
if (data != null) data = null;
if (imp != null) imp.close();
//数値の取得
getValue();
//キャンバスの元画像の取得と処理
imp0 = IJ.openImage("適当な画像ファイル");
size = imp0.getWidth();
//stack画像にする
for (int i = 0; i < number; i++) {
IJ.run(imp0, "Add Slice", "");
}
IJ.run(imp0, "Delete Slice", "");
imp = imp0.duplicate();
imp0.close();
ip = imp.getProcessor();
//データ類を入れる容器の作成
makeContainer();
}
////////////////////////////////////////////////
//JTextFieldからの数値の取得
public void getValue() {
dc = Float.parseFloat(dcField.getText());
cellSize = Integer.parseInt(cellField.getText());
obj = Integer.parseInt(objField.getText());
rate = Integer.parseInt(rateField.getText());
number = Integer.parseInt(timeField.getText());
}
////////////////////////////////////////////////
//容器の作成
public void makeContainer () {
//リセット
if (rawData != null) rawData = null;
if (data != null) data = null;
//初期化
data = new float [number][size][size];
localNumber = number*rate*2;
rawData = new float [localNumber][size][size];
//初期データの入力 (中心部以外は"0"に)
for (int i = 0; i < localNumber; i++) {
for (int j = 0; j < size; j++) {
for (int k = 0; k < size; k++) {
rawData[i][j][k] = 0.0f;
}
}
}
//中心部に値の導入
rawData [0][size/2][size/2] = 10000.0f;
rawData [0][size/2-1][size/2-1] = 10000.0f;
rawData [0][size/2-1][size/2] = 10000.0f;
rawData [0][size/2-1][size/2+1] = 10000.0f;
rawData [0][size/2][size/2-1] = 10000.0f;
rawData [0][size/2][size/2-1] = 10000.0f;
rawData [0][size/2+1][size/2-1] = 10000.0f;
rawData [0][size/2+1][size/2] = 10000.0f;
rawData [0][size/2+1][size/2+1] = 10000.0f;
}
////////////////////////////////////////////////
//各種計算
public void calc () {
//pixelサイズの計算
pixelSize = cellSize*1.0f/obj;
//伝播速度の計算
speed = dc*dt;
}
////////////////////////////////////////////////
//シミュレーションの計算
public void makeSimulate () {
//各種計算
calc();
//シミュレーションの計算
for (int i = 1; i < localNumber ; i++) {
for (int j = 1; j < size - 2; j++) {
for (int k = 1; k < size - 2; k++) {
rawData [i][j][k] = rawData[i-1][j][k] +
speed*(rawData[i-1][j-1][k] + rawData[i-1][j][k-1]+rawData[i-1][j][k+1] + rawData[i-1][j+1][k]
- 4*rawData[i-1][j][k])/(pixelSize*pixelSize);
}
}
}
//シミュレーション結果の記入 (イメージの撮影間隔を反映させる)
for (int i = 0; i < number; i++) {
for (int j = 0; j < size; j++) {
for (int k = 0; k < size; k++) {
data[i][j][k] = rawData[i*rate*2][j][k];
}
}
}
}
////////////////////////////////////////////////
//結果の表示
public void display() {
//シミュレーション画像の作成
for (int i = 0; i < number; i++) {
imp.setSlice(i+1);
for (int j = 0; j < size; j++) {
for (int k = 0; k < size; k++) {
ip.putPixelValue(j, k, data[i][j][k]);
}
}
}
//結果を見やすくする
IJ.run(imp, "Fire", "");
IJ.run("Brightness/Contrast...");
IJ.run(imp, "Enhance Contrast", "saturated=0.35");
//TimeとScale Barの表示
ip.setColor(Color.red);
Font font = new Font("Arial", Font.PLAIN, 20);
ip.setFont(font);
for (int i = 1; i < number + 1; i++) {
imp.setSlice(i);
int localTime = (int)Math.round((i-1)*rate);
if (localTime/1000 == 0) {
ip.drawString (localTime + " msec", 5, size);
} else if (localTime/1000 > 0) {
ip.drawString(localTime/1000 + " sec " + localTime %1000 + " msec", 5, size);
}
ip.drawString("20μm", (int)(size - 20/pixelSize), 25);
ip.drawLine((int)(size - 20/pixelSize - 10), 25, size -10, 25);
}
imp.show();
}
////////////////////////////////////////////////
//Close
public void closeAll() {
imp.close();
frame.dispose();
}
////////////////////////////////////////////////
//コントローラーの作成
public void run (String arg) {
Dimension screenSize = Toolkit.getDefaultToolkit().getScreenSize();
int sw = screenSize.width;
frame = new JFrame ("Fick Diffution");
frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
frame.setBounds (sw - 320, 50, 300, 320);
JLabel dcLabel = new JLabel ("<html>DC (μm<sup>2</sup>/s) = ");
JLabel cellLabel = new JLabel ("Cell Size (μm) =");
JLabel objLabel = new JLabel("Objective lens =");
JLabel rateLabel = new JLabel ("Frame rate (ms) = ");
JLabel timeLabel = new JLabel ("Time (frame) =");
JButton doneButton = new JButton ("Done");
JButton closeButton = new JButton ("Close");
dcField = new JTextField("20");
cellField = new JTextField ("16");
objField = new JTextField("60");
rateField = new JTextField ("30");
timeField = new JTextField("100");
doneButton.addActionListener(this);
closeButton.addActionListener(this);
doneButton.setActionCommand("Done");
closeButton.setActionCommand("Close");
dcLabel.setBounds(95, 10, 100, 30);
cellLabel.setBounds(80, 50, 100, 30);
objLabel.setBounds(75, 90, 150, 30);
rateLabel.setBounds(70, 130, 150, 30);
timeLabel.setBounds(90, 170, 100, 30);
dcField.setBounds(185, 12, 50, 30);
cellField.setBounds(185, 50, 50, 30);
objField.setBounds(185, 90, 50, 30);
rateField.setBounds (185, 130, 50, 30);
timeField.setBounds(185, 170, 50, 30);
doneButton.setBounds(110, 210, 80, 30);
closeButton.setBounds(110, 240, 80, 30);
JPanel pane = new JPanel ();
pane.setLayout(null);
pane.add(dcLabel);
pane.add(rateLabel);
pane.add(cellLabel);
pane.add(objLabel);
pane.add(dcField);
pane.add(rateField);
pane.add(cellField);
pane.add(objField);
pane.add(timeLabel);
pane.add(timeField);
pane.add(doneButton);
pane.add(closeButton);
frame.getContentPane().add(pane, BorderLayout.CENTER);
frame.setVisible(true);
}
////////////////////////////////////////////////
//ボタンを押した時の反応
public void actionPerformed(ActionEvent e) {
String cmd = e.getActionCommand();
if (cmd.equals("Done")) {
makeCanvas();
makeSimulate();
display();
} else if (cmd.equals("Close")) {
closeAll();
}
}
}このコードを実行すると、下記の様なコントローラーが立ち上がり、
Doneボタンを押すと、設定した条件での２次元拡散方程式のシミュレーションを行い、下の様に結果をstack画像として出力してくれる。
(これはビデオレートで観察した場合に見られる(はず)のの拡散の一コマである。)
少しコードの中身を見ていく。
(シミュレーション用のキャンバスの作成等の説明は
EclipseでImageJのPlugin作成 -拡散シミュレーション(ランダムウォーク)- - 生物屋さんのためのゼロからのプログラミングを参照)
このシミュレーションは顕微鏡での観察結果を想定しているので、コントローラー内の
DC = Cell Size = Objective lens = Frame rate = Time (frame) =
はイメージング条件の設定項目である。
"DC"は見たい分子の拡散係数(ここでは細胞内でのの拡散係数を使用)、"Cell Size"はカメラの画素サイズ、"Objective lens"は対物レンズの倍率、"Frame rate"は画像の撮像間隔、"Time"はイメージング時間を意味する。つまり、上の設定条件で「画素サイズ 16μm x 16μmのセンサーを持つカメラと60x対物レンズを使って、拡散係数が"Dc"である分子の動きを exposure time = 30 msecで100 frame分、Stream Acquisitionした」ことを意味する。
2次元拡散方程式は、
ここで、u(x, y, t)は時間tでの(x, y)座標での濃度、"Dc"は拡散係数、"Δt"はData間の時間間隔を示し、"Δx"と"Δy"は隣のピクセルとの距離を示すので、ここでは"Δx = Δy = pixelSize"とした。
また、拡散方程式の安定性条件は、
なので、16μm x 16μmのセンサーを持つカメラと60x対物レンズを使った場合に、この条件を満たす様に"Δt (コード上では"dt") = 0.0005" (0.5 msec)とした。（この場合の、上記の値は"0.14"である。）
そして、これらを踏まえて、
//pixelサイズの計算 pixelSize = cellSize*1.0f/obj; //伝播速度の計算 speed = dc*dt;
//シミュレーションの計算
for (int i = 1; i < localNumber ; i++) {
for (int j = 1; j < size - 2; j++) {
for (int k = 1; k < size - 2; k++) {
rawData [i][j][k] = rawData[i-1][j][k] +
speed*(rawData[i-1][j-1][k] + rawData[i-1][j][k-1]+rawData[i-1][j][k+1] + rawData[i-1][j+1][k]- 4*rawData[i-1][j][k])/(pixelSize*pixelSize);
}
}
}上記の部分でシミュレーションを行っている。
ここで得られた rawDataは、"0.5msec"毎の数値変化を収納していることになる。
("dt=0.5msec"なので)
そのため、
//シミュレーション結果の記入 (イメージの撮影間隔を反映させる)
for (int i = 0; i < number; i++) {
for (int j = 0; j < size; j++) {
for (int k = 0; k < size; k++) {
data[i][j][k] = rawData[i*rate*2][j][k];
}
}
}この部分で、コントローラー上のFrame rate (ms) を反映したデータを作っている。
(注："i*rate*2"の部分で時間の補正をした)
上で示したシミュレーション結果の画像はビデオレートのものだが、Frame rateを"100"に設定すると、
この様に、Frame間隔が長い場合のシミュレーションもすることが可能である。
また、初期値の入力時に、下記の様に
//初期データの入力
Random random = new Random();
for (int i = 0; i < localNumber; i++) {
for (int j = 0; j < size; j++) {
for (int k = 0; k < size; k++) {
if (i < 50 && random.nextInt(2500) == 0) {
rawData[i][j][k] = 10000.0f;
} else {
rawData[i][j][k] = 0;
}
}
}
}Random メソッドを使って、値が入る座標やframeをランダムになる様にすれば、
上の様に、色々な場所から様々なタイミングで拡散が始まる様子を観察できる。
