PythonでDICOM画像をなんとかするの続きです。
2017年09月時点での情報になります。
PyDicomは、現在githubからダウンロードすると、Versionが1.0.0a1になっているようです。
"""Pure python package for DICOM medical file reading and writing."""
__version__ = '1.0.0a1'
__version_info__ = (1, 0, 0, 'alpha', 0)
ただし、conda-forgeやPyPIでインストールできるpydicomは、まだ0.9.9になっているようです。1.0以降はいろいろと変更点があるので、ドキュメントにもわざわざTransition to pydicom 1.xという章が作られていますが、ひとまず、0.9.9準拠でこのページを作成することにします。
DICOM画像のサンプルを取得する
先に書いたPythonでDICOM画像をなんとかするでは、DICOM画像のサンプルの入手先としてOsirixのDICOM Image Libraryを紹介しましたが、現在ではOsirixユーザーでなければアクセスできなくなっているようです。
というワケで、DICOM画像のサンプルの入手先としては、以下の2つを紹介しておくことにします。
DICOMヘッダを全て表示する
知っているDICOMヘッダについては、以下のように表示することができます。これは、何のModalityで撮影された画像であるか？を表示します。
import dicom
dcm = dicom.read_file('CT_small.dcm')
print(dcm.Modality)
また、とにかく全部テキストで表示するだけであれば、特に工夫はいりません。
import dicom
dcm = dicom.read_file('CT_small.dcm')
print(dcm)
では、この print(dcm) のように、全てのヘッダを表示しよう(ただし、画像は表示しない)と考えたら、どうしたらよいのか？ということで、少しソースコードを追いかけてみました。
すると、pydicom/dataelem.pyにヒントがありました。
def __str__(self):
"""Return str representation of the element."""
repVal = self.repval
if self.showVR:
s = "%s %-*s %s: %s" % (str(self.tag), self.descripWidth,
self.description()[:self.descripWidth],
self.VR, repVal)
else:
s = "%s %-*s %s" % (str(self.tag), self.descripWidth,
self.description()[:self.descripWidth], repVal)
return s
どうやら、showVRというプロパティがDICOMヘッダのデータ型を表示させるかどうか？のフラグになっているようです。
VRには、どんな種類があるのか？というのは、DICOMの仕様を見るのが確実でしょう。
というところを踏まえて、dataelem.pyの__str__()をそのままパクると、全てのメタデータのタグと名前と型と表示をlistに.append()するなどして、全データを取得することができそうです。
import dicom
dcm = dicom.read_file('CT_small.dcm')
dcm_header_list = []
for data_element in dcm:
dcm_header = '{0}\t{1}\t{2}\n'.format(str(data_element.tag),
data_element.name,
data_element.repval)
dcm_header_list.append(dcm_header)
3D画像データからMIP画像を作成する
PythonでDICOM画像をなんとかするに書いたように、sliceを組み合わせて、3Dデータを作成することができます。
DICOMのファイルから画像データは、pixel_arrayからnumpyのarrayとして取り出すことができます。
import os
import dicom
import numpy as np
dcm_list = []
for root_d, _, files in os.walk('BRAINIX'):
for f in files:
_, ext = os.path.splitext(f)
if ext == '.dcm':
dcm_list.append(os.path.join(root_d, f))
ref_dcm = dicom.read_file(dcm_list[0])
width = ref_dcm.Columns
height = ref_dcm.Rows
depth = len(ref_dcm)
img3d = np.zeros((depth, height, width),
dtype = ref_dcm.pixel_array.dtype)
for dcm_file in dcm_list:
idx = dcm_list.index(dcm_file)
dcm = dicom.read_file(dcm_file)
img3d[idx, :, :] = dcm.pixel_array()
これにより、img3dという三次元のアレイにDICOMの画像データが格納されます。
これを、可視化するために、MIP(Maximum Intensity Projection)あるいはMinIP(Minimum Intensity Projection)という方法を使って、二次元の画像に変換します。
三次元の画像データをどの軸に沿った断面の画像を作るのか？、そして、その軸に沿った二次元画像を作成する際に、どの軸の画素を代表値として用いるのか？ということですが、MIPは最も輝度が高い画素を、MinIPは逆に最も輝度が低い画素を代表値として用います。
また、人体に対してどの断面を使うのか？ということについては、それぞれの断面に名前がついています。Wikipediaの記事: Anatomical_plane
最大値を使うにせよ、最小値を使うにせよ、ここからはnumpy arrayの使いこなしになります。
# 前略
max0 = np.nanmax(img3d, axis=0)
max1 = np.nanmax(img3d, axis=1)
max2 = np.nanmax(img3d, axis=2)
今回のBRAINIXのMRI画像の場合は、axis=0がCoronal Planeにあたるので、max0には、Coronal PlaneのMIP画像が格納されます。あとは、matplotlib.imshow()なり、 PILでim = Image.fromarray(max0)してからim.show()するなりすればOKだと思います。
元の画像の断面がaxialなのかcoronalなのかsagittalなのかは、元の画像から判断するしかなさそうだけれど、こうしてimg3dのような三次元アレイにしてしまえば、どの断面から切り出すことも可能になります。
一般論としては、CTはaxial planeで撮像するので、[axial, coronal, sagittal]の三次元になっているように思います。
まとめ
PythonでDICOM画像をなんとかするの続きを書いてみました。