X線投影シミュレーションソフトウェア ～gVirtualXrayの使い方１～

Skull-san

今回はgVirtualXrayというXpやCTなどが可能なシミュレーションソフトを紹介します。

gVirtualXRay (gVXR)は、GPU上で動作するオープンソースのX線イメージシミュレーションフレームワークです。X線の減弱をビール–ランバートの法則（I＝Ioe^（-μx）の式）に基づいてモデル化し、3Dオブジェクト（ポリゴンメッシュ）を通過する光子の吸収を計算します​

gVXRは減弱現象のみをモデル化した解析的シミュレーションです（散乱線やノイズの効果は含まない）。

そのため散乱の影響が小さい条件においては、非常に現実的な画像を高速に得ることができます。必要に応じて散乱や量子ノイズの効果を後処理で付加することも可能です。

一方、gVXRには3Dシーンのリアルな表示機能（疑似X線画像だけでなく、X線源・検出器配置や物体の3次元位置関係を可視化する機能）も備わっており、OpenGLを用いたインタラクティブな3Dビューアでシミュレーション環境を確認できます​。

単純X線撮影画像だけでなく、被写体を多角度から透視した投影像を連続的に生成することでCTスキャンのシミュレーションも可能です。

また、X線源のエネルギースペクトルを単色（単一エネルギー）だけでなく任意の分布で設定でき、スペクトルを考慮したシミュレーション（単色X線CTやスペクトラルCT）も行えます。

被写体となる3Dオブジェクトはポリゴンメッシュで表現し、それぞれに材質（物質情報）を割り当てることができます。材質指定は柔軟で、元素記号や化合物名による定義のほか、複数の元素からなる混合物、密度、CT値（HU）で設定することも可能です。

Skull-san

本記事で紹介するコードはこのページのFirst simulation X-ray imageのipynbをもとに作成しています。
gvxrはインタラクティブな関数を提供しているはずですが、上手く動作しなかったので（よくわからなかったので）、少し強引な方法でやっています。
なお、環境はWin10です。

Google colabratoryの準備

Google colab.を使うことによってほとんど準備を必要とせずにpythonを使える環境が手に入ります。

ネットに膨大な情報がありますので、ここでgoogle colab.の使い方を紹介する必要はないでしょう。このサイトを見て、準備してください。

Googleアカウントさえあればすぐに使えます。

インストール

Google colab.で以下のコマンドを打ちます。

!pip install matplotlib numpy ipywidgets gvxr
!jupyter nbextension enable --py widgetsnbextension

次に必要なライブラリをimportします。

import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import PowerNorm
from ipywidgets import interact, FloatSlider
import ipywidgets as widgets
from IPython.display import display
from IPython.display import clear_output
from gvxrPython3 import gvxr

大本のGVXR（コア部分）とSimple GVXR（ラッパー）のバージョンが違うとトラブルが起こるかも知れないので、一応printして確認します。一致してればOKです。

print(gvxr.getVersionOfCoreGVXR())  # Core GVXRのバージョン
print(gvxr.getVersionOfSimpleGVXR())  # Simple GVXRのバージョン

次にOpenGLのcontextを作成します。まぁ準備ですので、気にしないでください。

gvxr.createOpenGLContext()

画像を出力したりする場合に備えてoutput用のフォルダを作成します。
※本記事では出力していないので、ここは無視してもらっていいです。自分用のメモ書きです。

# 出力フォルダの設定
output_path = "./output"
os.makedirs(output_path, exist_ok=True) #存在していればスキップ、存在しなければ新しく作成する

投影する用のファントム（STLファイルをアップロードしておきます。）
何でもいいです。私は手元にあった人体骨3DのSTLファイルをアップロードしました。D＆DでOK。

では、準備が整ったので処理の中身に入っていきます。

# 初期設定
def initialize_scene():
    gvxr.setSourcePosition(-40.0,  0.0, 0.0, "cm")
    gvxr.usePointSource()
    gvxr.setMonoChromatic(80, "keV", 1000)
    gvxr.setDetectorPosition(10.0, 0.0, 0.0, "cm")
    gvxr.setDetectorUpVector(0, 0, -1)
    gvxr.setDetectorNumberOfPixels(320, 900)
    gvxr.setDetectorPixelSize(0.5, 0.5, "mm")

    # STLファイル読み込み
    fname = "./human2.stl"
    if not os.path.exists(fname):
        raise IOError(fname)

    gvxr.loadMeshFile("Human", fname, "mm")
    gvxr.moveToCentre("Human")

    gvxr.setCompound("Human", "H2O")
    gvxr.setDensity("Human", 1.0, "g/cm3")

initialize_scene()

コードを見て分かるように、
線源の位置、エネルギー、検出器の位置・向き、ピクセル数・サイズを指定可能です。上記の設定は初期設定でこの状態からインタラクティブに変更していきます。

また、STLファイルの物質は水として扱い、減弱を計算します。
initialize_scene()という関数を定義して、呼んでいます。

次に、画像をアップデートするための関数を作成していきます。おそらくこんな方法を取らなくてもいいはずなんだが・・・。

def update_simulation(src_x, src_y, src_z, det_x, det_y, det_z, up_x, up_y, up_z, gamma):
    # 出力をクリア（毎回表示を更新）
    clear_output(wait=True)

    # 線源位置を更新
    gvxr.setSourcePosition(src_x, src_y, src_z, "cm")
    gvxr.setDetectorPosition(det_x, det_y, det_z, "cm")
    gvxr.setDetectorUpVector(up_x, up_y, up_z)

    # X線画像を再計算
    x_ray_image = gvxr.computeXRayImage()

    # 3Dシーンを更新
    gvxr.displayScene()

    # 画像表示
    plt.imshow(np.array(x_ray_image), cmap="gray", norm=PowerNorm(gamma=gamma))
    plt.title("X-ray Image")
    plt.colorbar()
    plt.axis("off")
    plt.show()

    # 再度UI表示（消えないように）
    display(ui)

最後にパラメータを調整するスライダーを作成して、画像化までのコードです。

# スライダー群を定義
src_x = widgets.FloatSlider(min=-100, max=0, step=1, value=-40, description="Source X (cm)")
src_y = widgets.FloatSlider(min=-50, max=50, step=1, value=0, description="Source Y (cm)")
src_z = widgets.FloatSlider(min=-50, max=50, step=1, value=0, description="Source Z (cm)")

det_x = widgets.FloatSlider(min=0, max=100, step=1, value=10, description="Detector X (cm)")
det_y = widgets.FloatSlider(min=-50, max=50, step=1, value=0, description="Detector Y (cm)")
det_z = widgets.FloatSlider(min=-50, max=50, step=1, value=0, description="Detector Z (cm)")

up_x = widgets.FloatSlider(min=-1, max=1, step=0.1, value=0, description="Up Vector X")
up_y = widgets.FloatSlider(min=-1, max=1, step=0.1, value=0, description="Up Vector Y")
up_z = widgets.FloatSlider(min=-1, max=1, step=0.1, value=-1, description="Up Vector Z")

gamma = widgets.FloatSlider(min=0.1, max=3.0, step=0.1, value=1.0, description="Gamma")

# 実行ボタン
run_button = widgets.Button(description="投影", button_style="success")


# 実行ボタンが押された時の処理
def on_run_button_clicked(b):
    update_simulation(
        src_x.value, src_y.value, src_z.value,
        det_x.value, det_y.value, det_z.value,
        up_x.value, up_y.value, up_z.value,
        gamma.value
    )

# ボタンにイベントを関連づけ
run_button.on_click(on_run_button_clicked)

# スライダーとボタンを表示
ui = widgets.VBox([
    src_x, src_y, src_z,
    det_x, det_y, det_z,
    up_x, up_y, up_z,
    gamma,
    run_button
])
display(ui)

これを実行すると以下のようなスライダーが表示されますので、調節して「投影」を押します。

以下の画像はいくつかの条件における、投影画像です。

jupyter notebookだと以下のように表示されます。

線源の距離を近くすると拡大され、線源・検出器の方向を変えると画が流れているのが分かります。

必要ないらしいですが、一応、後片付けをしてプログラムを終了します。

# シミュレーション終了時にコンテキストをクリーンアップ
gvxr.destroy()
print("Simulation ended!")

Skull-san

今回はここまでとします。
これを利用して、自作のソフトウェアを作ることも可能です。