本記事には広告が含まれています。
This article contains advertisements.

Python

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

 

Skull-san

前回紹介したgVirtualXrayを使ってCTのシミュレーションをやっていきます。再構成にはCIL（Core Imaging Library）を使用します。ただし、google colab.やWindowsではCILのインストールがうまくいかなかったので今回はWin10上で動作するWSL（Ubuntu）を使用していきます。GPUはNDIVIA GeForce 1060を使っています。

私が行った流れとしては・・・
１．WSLの準備
２．WSLでGPUを使うためにCUDAドライバをインストール
３．VSCodeからWSLを操作するための準備
４．Minicondaをインストール
５．CILのインストール、仮想環境（venv_CIL）を作成
６．５で作成した仮想環境をactivate
７．gVirtualXrayのインストール
８．3D表示用にK3Dをインストール（結局使っていないけど）
９．作成したPythonコードの実行
という感じです。

Linux環境を用意するにあたってvirtualBoxなどではGPUは容易に使えません（方法はあるようです）。また、VSCodeも必須ではありませんので不要であれば飛ばしてください。

Skull-san

本記事で紹介するコードはこの論文に書いてあるコードをもとに作成しています。gvxrは高速ですが、CILの再構成（特にSIRT）はそこそこ時間がかかります。

もしかしたら色々と試行錯誤している途中で追加インストールしたものがあるかも知れません。記事では書き忘れているものがあるかもしれませんので、エラーが出たらその都度対応してください。


事前準備

１～４（WSL、VSCode、Miniconda）の準備は先ほどのリンクを参照しながら準備してください。


ここまで終わっていることを想定して、次に進みます。

gvxrとCILのインストール

上記の準備が終わって、VSCodeを実行した様子です。
WSLに接続しており、(base)と表示されています。Minicondaをインストールしたらデフォルトのbase環境が立ち上がっています。

開発したいものに応じて、色々なパッケージをインストールしますが、普通はこのbaseにはインストールせずに、仮想環境（venvなどと名前をつける）をそのつど作成して、そこに環境構築をしていきます。

では、まずCILのgithubにあるコマンドをターミナルに打ってみましょう。

conda create --name venv_CIL -c conda-forge -c https://software.repos.intel.com/python/conda -c ccpi cil=24.3.0 ipp=2021.12 astra-toolbox=*=cuda* tigre ccpi-regulariser tomophantom ipykernel ipywidgets scikit-image


環境が作成されたら、指示されたコマンドを打ってactivateします。

今回は仮想環境の名前がvenv_CILなので(–nameの後の文字)、以下のようになります。

conda activate venv_CIL

次にgVirtualXrayをgitからcloneします。
Condaなどの環境変数の影響を完全に無視してGitを実行するためのコマンドも含んでいます。（こうする必要があったのは、私が色々試していたせいかもしれません・・・・）

env -i HOME=$HOME PATH=/usr/bin:/bin git clone --recursive https://github.com/effepivi/gVXR2CIL.git

無事にできたら、以下のようにディレクトリ（gVXR2CIL）ができています。

以下のようにまだgVirtualXrayは使えません。

なので

pip install gvxr

を打ってからimportします。

正常にimportできました。これで準備環境です。

使用する３Dデータ

ここから貼るコードをすべてつなげて、実行すれば最後まで進むはずです（というか、ほとんどこの論文のコードです）。本記事の最後に、全コードを載せますのでそちらを見ていただければOKです。なお、ChatGPTにコードの整理およびコメントによる説明を追加してもらいました。

流れとしては前半がGVXRによるprojection画像の取得、後半がCILによる画像再構成です。再構成法はFDKとSIRTです。FDKはfeldkampらの解析的な方法、SIRTは逐次近似法です。

今回も前回同様に人の骨の3Dモデル（human2.stl）を使用します。作業ディレクトリにいれておきます。読者の皆さんは各自、好きな3Dデータを使用してください。

こんなデータです。

gVirtualXrayによるprojection画像の取得

まず最初に必要なものをインポートします。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import os
import numpy as np
import matplotlib.pyplot as plt
from gvxrPython3 import gvxr
from cil.io import TIFFStackReader, TIFFWriter
from cil.framework import AcquisitionGeometry, AcquisitionData
from cil.recon import FDK
from cil.optimisation.algorithms import SIRT
from cil.optimisation.functions import IndicatorBox
from cil.plugins.astra.operators import ProjectionOperator

ここでnumpy等、エラーが出る場合にはインストールしておいてください。

# -------------------------------
# 1. gVXRによるCTスキャンシミュレーション
# -------------------------------
print("OpenGLウィンドウを作成中...")
gvxr.createWindow()
gvxr.setWindowSize(500, 500)

# X線源と検出器の設定
gvxr.setSourcePosition(-100.0, 0.0, 0.0, "cm")
gvxr.usePointSource()
gvxr.setMonoChromatic(60.0, "keV", 1000)
gvxr.setDetectorPosition(100.0, 0.0, 0.0, "cm")
gvxr.setDetectorUpVector(0, 0, 1)
gvxr.setDetectorNumberOfPixels(512, 512)
gvxr.setDetectorPixelSize(0.5, 0.5, "mm")
gvxr.setScintillator("CsI", 0.5, "mm")

# 対象物（human2.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")

# 投影取得・保存
gvxr.setZoom(2)
gvxr.displayScene()
screenshot = gvxr.takeScreenshot()

output_dir = "ct_projections2"
os.makedirs(output_dir, exist_ok=True)
number_of_projections = 180

gvxr.computeCTAcquisition(
    os.path.join(output_dir, f"projections-{number_of_projections}"),
    os.path.join(output_dir, f"screenshots-{number_of_projections}"),
    number_of_projections, 0.0, True, 360.0, 60,
    0.0, 0.0, 0.0, "mm",
    *gvxr.getDetectorUpVector(), True, 0
)

print("CTスキャンシミュレーション完了。")

コードを読んだら分かりますが、
X線管球はX軸のマイナス方向、検出器はX軸のプラス方向に配置しています。
なので、マイナス方向からプラス方向に向けて投影されることになります。

エネルギーは単一の60 keVです。投影方向は180なので、180枚の画像が出力されます。

projectionデータはカレントディレクトリのct_projections2に出力されます。ディレクトリが無い場合は作られます。

このコードを実行した結果、得られる画像は以下のようなものです。

CILによる再構成

では、ここから後半になります。

# -------------------------------
# 2. 投影画像の読込と吸収変換
# -------------------------------
print("TIFFStackReaderで投影データを読み込み中...")
reader = TIFFStackReader(file_name=os.path.join(output_dir, f"projections-{number_of_projections}"), dtype=np.float32)
data_original = reader.read()

print("投影データを正規化・ログ変換します...")
data_normalised = data_original / data_original.max()
data_normalised[data_normalised < 1e-9] = 1e-9

data_absorption = -np.log(data_normalised)

ここでは、先ほど出力したprojection画像（TIFF形式）を読み込んでいます。

また、正規化およびログ変換を（非ゼロ化も）しています。

CTの教科書に（非常に重要！！！と）書いてあるようにX線の物質中の吸収はI=I0exp(－μx）で表されます。既知の情報はIとI0だけであり、式を変形(μx = -ln(I/I0))して吸収量（μx）を求めます。よって対数変換という処理が必要です。

次のコードは各投影が、どこからどの方向に撮られたものかをCILに伝える部分になります。

# -------------------------------
# 3. AcquisitionGeometryの作成
# -------------------------------
print("AcquisitionGeometryを構築中...")
source_position_cm = gvxr.getSourcePosition("cm")
detector_position_cm = gvxr.getDetectorPosition("cm")
detector_dir_x = gvxr.getDetectorRightVector()
detector_dir_y = gvxr.getDetectorUpVector()
rotation_axis_pos = (0, 0, 0)
rotation_axis_dir = gvxr.getDetectorUpVector()
angles_degs = np.linspace(0, 360, number_of_projections, endpoint=True)

geometry = AcquisitionGeometry.create_Cone3D(
    source_position=source_position_cm,
    detector_position=detector_position_cm,
    detector_direction_x=detector_dir_x,
    detector_direction_y=detector_dir_y,
    rotation_axis_position=rotation_axis_pos,
    rotation_axis_direction=rotation_axis_dir
)
geometry.set_angles(angles_degs, angle_unit="degree")
geometry.set_panel(
    gvxr.getDetectorNumberOfPixels(),
    gvxr.getDetectorPixelSpacing("cm")
)
geometry.set_labels(["angle", "vertical", "horizontal"])

# gVXRのウィンドウを終了
gvxr.terminate()

次からいよいよ再構成です。

# -------------------------------
# 4. CILによる再構成処理
# -------------------------------
# AcquisitionData作成・並べ替え
data = AcquisitionData(data_absorption, deep_copy=False, geometry=geometry)
data.reorder(order="tigre")

# FDK再構成
print("FDK再構成を実行します...")
image_geometry = data.geometry.get_ImageGeometry()
fdk = FDK(data, image_geometry)
recon_fdk = fdk.run()

# TIFF保存
fdk_output_dir = os.path.join(output_dir, "recon-FDK")
os.makedirs(fdk_output_dir, exist_ok=True)
TIFFWriter(data=recon_fdk, file_name=os.path.join(fdk_output_dir, "recon_fdk")).write()
print("FDK再構成結果を保存しました。")

# SIRT再構成
print("SIRT再構成を実行します...")
data.reorder(order="astra")
A = ProjectionOperator(image_geometry, geometry, "gpu")
x0 = image_geometry.allocate()
constraint = IndicatorBox(lower=0)

sirt = SIRT(initial=x0, operator=A, data=data, constraint=constraint, max_iteration=500)
sirt.update_objective_interval = 50
sirt.run(500)
recon_sirt = sirt.solution

# TIFF保存
sirt_output_dir = os.path.join(output_dir, "recon-SIRT")
os.makedirs(sirt_output_dir, exist_ok=True)
TIFFWriter(data=recon_sirt, file_name=os.path.join(sirt_output_dir, "recon_sirt")).write()
print("SIRT再構成結果を保存しました。")

FDKはそんなに時間かかりませんが、SIRT（500 iteration）は結構時間がかかります。

再構成した画像を保存している場所に行くにはexplorerに

\\wsl$\Ubuntu\home\test

を貼り付けます。

こんなかんじです。

開くと再構成画像が保存されているのが分かると思います。

では、再構成画像を見てみましょう。

FDK（128×128に解像度を下げています）

次がSIRT

FDKよりもストリークアーチファクトが無くてきれいですね。

コード全文

最後に全コードを貼って終わりになります。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
CTシミュレーション＋FDK/SIRT再構成コード（統合版）

【処理の流れ】
1. STLファイル（human2.stl）を読み込み、gVirtualXray (gvxr) で投影画像をシミュレーション
2. 投影画像をTIFFとして保存
3. CILで投影画像を読み込み、透過画像を吸収画像へ変換
4. CT幾何情報を構築（source/detector位置、角度、画素サイズなど）
5. FDK再構成を実行、TIFFに保存
6. SIRT再構成を実行、TIFFに保存
"""

import os
import numpy as np
import matplotlib.pyplot as plt
from gvxrPython3 import gvxr
from cil.io import TIFFStackReader, TIFFWriter
from cil.framework import AcquisitionGeometry, AcquisitionData
from cil.recon import FDK
from cil.optimisation.algorithms import SIRT
from cil.optimisation.functions import IndicatorBox
from cil.plugins.astra.operators import ProjectionOperator

# -------------------------------
# 1. gVXRによるCTスキャンシミュレーション
# -------------------------------
print("OpenGLウィンドウを作成中...")
gvxr.createWindow()
gvxr.setWindowSize(500, 500)

# X線源と検出器の設定
gvxr.setSourcePosition(-100.0, 0.0, 0.0, "cm")
gvxr.usePointSource()
gvxr.setMonoChromatic(60.0, "keV", 1000)
gvxr.setDetectorPosition(100.0, 0.0, 0.0, "cm")
gvxr.setDetectorUpVector(0, 0, 1)
gvxr.setDetectorNumberOfPixels(512, 512)
gvxr.setDetectorPixelSize(0.5, 0.5, "mm")
gvxr.setScintillator("CsI", 0.5, "mm")

# 対象物（human2.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")

# 投影取得・保存
gvxr.setZoom(2)
gvxr.displayScene()
screenshot = gvxr.takeScreenshot()

output_dir = "ct_projections2"
os.makedirs(output_dir, exist_ok=True)
number_of_projections = 180

gvxr.computeCTAcquisition(
    os.path.join(output_dir, f"projections-{number_of_projections}"),
    os.path.join(output_dir, f"screenshots-{number_of_projections}"),
    number_of_projections, 0.0, True, 360.0, 60,
    0.0, 0.0, 0.0, "mm",
    *gvxr.getDetectorUpVector(), True, 0
)

print("CTスキャンシミュレーション完了。")

# -------------------------------
# 2. 投影画像の読込と吸収変換
# -------------------------------
print("TIFFStackReaderで投影データを読み込み中...")
reader = TIFFStackReader(file_name=os.path.join(output_dir, f"projections-{number_of_projections}"), dtype=np.float32)
data_original = reader.read()

print("投影データを正規化・ログ変換します...")
data_normalised = data_original / data_original.max()
data_normalised[data_normalised < 1e-9] = 1e-9

data_absorption = -np.log(data_normalised)

# -------------------------------
# 3. AcquisitionGeometryの作成
# -------------------------------
print("AcquisitionGeometryを構築中...")
source_position_cm = gvxr.getSourcePosition("cm")
detector_position_cm = gvxr.getDetectorPosition("cm")
detector_dir_x = gvxr.getDetectorRightVector()
detector_dir_y = gvxr.getDetectorUpVector()
rotation_axis_pos = (0, 0, 0)
rotation_axis_dir = gvxr.getDetectorUpVector()
angles_degs = np.linspace(0, 360, number_of_projections, endpoint=True)

geometry = AcquisitionGeometry.create_Cone3D(
    source_position=source_position_cm,
    detector_position=detector_position_cm,
    detector_direction_x=detector_dir_x,
    detector_direction_y=detector_dir_y,
    rotation_axis_position=rotation_axis_pos,
    rotation_axis_direction=rotation_axis_dir
)
geometry.set_angles(angles_degs, angle_unit="degree")
geometry.set_panel(
    gvxr.getDetectorNumberOfPixels(),
    gvxr.getDetectorPixelSpacing("cm")
)
geometry.set_labels(["angle", "vertical", "horizontal"])

# gVXRのウィンドウを終了
gvxr.terminate()

# -------------------------------
# 4. CILによる再構成処理
# -------------------------------
# AcquisitionData作成・並べ替え
data = AcquisitionData(data_absorption, deep_copy=False, geometry=geometry)
data.reorder(order="tigre")

# FDK再構成
print("FDK再構成を実行します...")
image_geometry = data.geometry.get_ImageGeometry()
fdk = FDK(data, image_geometry)
recon_fdk = fdk.run()

# TIFF保存
fdk_output_dir = os.path.join(output_dir, "recon-FDK")
os.makedirs(fdk_output_dir, exist_ok=True)
TIFFWriter(data=recon_fdk, file_name=os.path.join(fdk_output_dir, "recon_fdk")).write()
print("FDK再構成結果を保存しました。")

# SIRT再構成
print("SIRT再構成を実行します...")
data.reorder(order="astra")
A = ProjectionOperator(image_geometry, geometry, "gpu")
x0 = image_geometry.allocate()
constraint = IndicatorBox(lower=0)

sirt = SIRT(initial=x0, operator=A, data=data, constraint=constraint, max_iteration=500)
sirt.update_objective_interval = 50
sirt.run(500)
recon_sirt = sirt.solution

# TIFF保存
sirt_output_dir = os.path.join(output_dir, "recon-SIRT")
os.makedirs(sirt_output_dir, exist_ok=True)
TIFFWriter(data=recon_sirt, file_name=os.path.join(sirt_output_dir, "recon_sirt")).write()
print("SIRT再構成結果を保存しました。")

print("すべての処理が完了しました。")


Skull-san

読んでいただきありがとうございました～。

RELATED POST

COMMENT

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です

CAPTCHA