pythonとTIGREを使ったコーンビームCT（CBCT）の再構成

#はじめに
色々ありまして3次元のCBCT(Cone-Beam Computed Tomography)を投影画像から再構成したくなりました。しかし日本語で「CBCT　再構成　python」と検索しても、再構成の理論についてはヒットするのですが、コードでどうやるかについてはあまり出てきません。特に2次元の再構成は多少なりとも見つかるのですが、3次元の再構成はほとんど見つかりません。
そこで今回、TIGRE（https://github.com/CERN/TIGRE ）というライブラリを使ってpythonで3次元CBCTを再構成しましたので、そのやり方を共有します。

####※注意事項
本記事はあくまでTIGREを使った3次元CBCTの再構成方法を解説するだけであって、再構成プログラムを0から作成するわけではありません。またCBCT再構成の理論的な解説はありません。
理論とかプログラムの詳細とかはどうでもよく、とにかく3次元CBCTを再構成できればわしゃそれでいいんじゃ！！という人向けです。ご了承ください。
ちなみにCBCTについて詳しく知りたい方にはコーンビームCT画像再構成の基礎¹がおすすめです。

#PC環境

• OS : Windows 10 Professional
• 言語 : Python 3.7.3
• GPU : Geforce RTX 2080 Ti
• Python : 3.7.3
• CUDA : 10.1
• MSVS : 2017
• TIGRE : 2.0

#CTとCBCTの再構成の違い
すごーーーくざっくり説明すると、CTとCBCTでは検出器に入ってくるX線がファンビームか、コーンビームかという点が大きく異なります。そのため再構成手法も変わってきます。
CTでは検出器に対してX線がほぼ垂直に入射してくるため（＝ファンビーム）、フィルター逆投影法（FBP法: Filtered Back Projection）という手法で再構成できます。一方でCBCTの場合は検出器に対してX線が角度（コーン角）を持って入射してくるため（＝コーンビーム）、FBP法ではなくFeldkampらが提案したFDK法²で再構成する必要があります。
細かい解説はコーンビームCT画像再構成の基礎¹を見てもらうとして、ここではとりあえず「CBCTを再構成する際は、FBP法ではなくFDK法を使う」ということだけ覚えていてください。

#TIGREについて
今回はTIGRE(https://github.com/CERN/TIGRE) というCERN（欧州原子核研究機構）の方たちが開発したライブラリを使って3次元CBCT再構成をしていきます。CERNの方は何でもできてすごいですね。
PythonのCBCT再構成のライブラリは検索すると他にもいくつか出てくるのですが、個人的にTIGREが一番簡単にインストールできて、かつGPUにも対応しているので使いやすいと思います。

#TIGREのインストール
Github上のInstallation Instructions for PythonのSimple Instructionsの手順通りに進めます。まずMSVC、CUDA、pythonをインストールした後、下記を実行します。

git clone https://github.com/CERN/TIGRE.git
#gitで取得したTIGRE/Python/に移動したあと
python setup.py install --user

setup.pyがエラーなく終了したら、TIGRE/Pythonフォルダにある、example.pyを実行してみてください。以下のような画像が出てきたら正しくインストールされています。

####example.pyの実行結果

####インストール時のエラー
ちなみに私の環境だとsetup.pyを実行した際に以下のようなエラーがでました。

nvcc fatal : Unsupported gpu architecture 'compute_86'

これは私の持っているRTX 2080 Tiがcompute_86アーキテクチャに対応していないのが原因？のようで、setup.pyの一部をコメントアウトすることで解決しました。

setup.py

#18行目あたり
COMPUTE_CAPABILITY_ARGS = [  # '-gencode=arch=compute_20,code=sm_20', #deprecated
    #'-gencode=arch=compute_30,code=sm_30',#deprecated
    '-gencode=arch=compute_37,code=sm_37',
    '-gencode=arch=compute_52,code=sm_52',
    '-gencode=arch=compute_60,code=sm_60',
    '-gencode=arch=compute_61,code=sm_61',
    '-gencode=arch=compute_70,code=sm_70',
    '-gencode=arch=compute_75,code=sm_75',
    #'-gencode=arch=compute_86,code=sm_86',#←これ
    '--ptxas-options=-v', '-c',
    '--default-stream=per-thread',
    ]

#CBCTの投影と再構成
TIGREをインストールできたらあとはもう簡単です。というのもexample.pyですでに投影画像の取得とFDK法による再構成ができています。ただしexample.pyではCBCTのGeometryがデフォルトの状態なので、Geometryを変更できるようにします。

FDK.py

from __future__ import print_function
from __future__ import division
import numpy as np
from matplotlib import pyplot as plt
import tigre
import tigre.algorithms as algs
from tigre.demos.Test_data import data_loader
from tigre.utilities.Measure_Quality import Measure_Quality

# Geometryの設定
geo = tigre.geometry(mode='cone',  default=False)
geo.DSD = 1500                                     # Distance Source Detector      (mm)
geo.DSO = 1000                                     # Distance Source Origin        (mm)
# Detector parameters
geo.nDetector = np.array((300, 300))               # number of pixels              (px)
geo.dDetector = np.array((0.8, 0.8))               # size of each pixel            (mm)
geo.sDetector = geo.nDetector * geo.dDetector    # total size of the detector    (mm)
# Image parameters
geo.nVoxel = np.array((256, 256, 256))             # number of voxels              (vx)
geo.sVoxel = np.array((128, 128, 128))             # total size of the image       (mm)
geo.dVoxel = geo.sVoxel/geo.nVoxel               # size of each voxel            (mm)
# Offsets
geo.offOrigin = np.array((0, 0, 0))                # Offset of image from origin   (mm)
geo.offDetector = np.array((0, 0))                 # Offset of Detector            (mm)

# Auxiliary
geo.accuracy = 0.5                                 # Accuracy of FWD proj          (vx/sample)
# Mode
geo.mode = 'cone'                                  # parallel, cone                ...
geo.filter = 'hann'                       #  None, shepp_logan, cosine, hamming, hann

#撮影角度, 0~360度での撮影
nangles = 360
angles = np.linspace(0, 2 * np.pi, nangles, endpoint=False, dtype=np.float32)

#デフォルトで用意されている頭部ファントムを取得
head = data_loader.load_head_phantom(geo.nVoxel)

#頭部ファントムの投影画像を取得
proj = tigre.Ax(head,geo,angles)

#FDKでの再構成
fdkout = algs.fdk(proj,geo,angles)

#描画
plt.subplot(1,3,1)
plt.imshow(fdkout[geo.nVoxel[0]//2])
plt.subplot(1,3,2)
plt.imshow(fdkout[:, geo.nVoxel[1]//2, :])
plt.subplot(1,3,3)
plt.imshow(fdkout[:, :, geo.nVoxel[2]//2])
plt.show()

#tigreで用意されている描画メソッド
tigre.plotProj(proj)
tigre.plotImg(fdkout)

####実行結果
FDK法による再構成画像です。

また、tigre.plotProj(proj)の部分で、角度毎の投影画像を確認できます。

####入力画像の変更
入力とするCT画像はnumpy.ndarray形式なので、別のCT画像で計算したい場合は入力画像をnumpy.ndarray形式に変更したのち、以下の部分を置き換えるだけです。

head = data_loader.load_head_phantom(geo.nVoxel)

また自分で用意した投影画像を使って再構成したい場合も同じく、投影画像をnumpy.ndarray形式に変えたのち、以下の部分を置き換えるだけです。

proj = tigre.Ax(head,geo,angles)

#まとめ
pythonで3次元CBCTの再構成をするため、TIGREというライブラリの使い方を解説しました。
CBCT再構成をやりたいけどよくわからん、という人の役に立てば幸いです。

※間違い等ありましたらご指摘お願いいたします。

####参考文献

1. コーンビーム再構成の基礎 (http://www.iryokagaku.co.jp/frame/03-honwosagasu/433/433.html) ↩ ↩²

2. https://doi.org/10.1364/JOSAA.1.000612 ↩