Python
dicom

PythonでDICOM画像をなんとかする その2

More than 1 year has passed since last update.

PythonでDICOM画像をなんとかするの続きです。
2017年09月時点での情報になります。

PyDicomは、現在githubからダウンロードすると、Versionが1.0.0a1になっているようです。

_versions.py
"""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にヒントがありました。

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画像をなんとかするの続きを書いてみました。