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

[PythonでDICOM画像をなんとかする](http://qiita.com/fukuit/items/ed163f9b566baf3a6c3f)の続きです。
2017年09月時点での情報になります。

[PyDicom](https://github.com/pydicom/pydicom)は、現在githubからダウンロードすると、Versionが1.0.0a1になっているようです。

```Python:_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](http://pydicom.readthedocs.io/en/stable/transition_to_pydicom1.html)という章が作られていますが、ひとまず、0.9.9準拠でこのページを作成することにします。

# DICOM画像のサンプルを取得する
先に書いた[PythonでDICOM画像をなんとかする](http://qiita.com/fukuit/items/ed163f9b566baf3a6c3f)では、DICOM画像のサンプルの入手先としてOsirixの[DICOM Image Library](http://www.osirix-viewer.com/resources/dicom-image-library/)を紹介しましたが、現在ではOsirixユーザーでなければアクセスできなくなっているようです。

というワケで、DICOM画像のサンプルの入手先としては、以下の2つを紹介しておくことにします。

+ [ImageJのExamples](https://imagej.nih.gov/ij/images/)
+ [PyDicomのtest_files](https://github.com/pydicom/pydicom/tree/master/pydicom/data/test_files)


# DICOMヘッダを全て表示する
知っているDICOMヘッダについては、以下のように表示することができます。これは、何のModalityで撮影された画像であるか？を表示します。

```Python
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にヒントがありました。

```Python: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の仕様](http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_6.2.html)を見るのが確実でしょう。

というところを踏まえて、```dataelem.py```の```__str__()```をそのままパクると、全てのメタデータのタグと名前と型と表示をlistに.append()するなどして、全データを取得することができそうです。

```Python
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画像をなんとかする](http://qiita.com/fukuit/items/ed163f9b566baf3a6c3f)に書いたように、sliceを組み合わせて、3Dデータを作成することができます。
DICOMのファイルから画像データは、pixel_arrayからnumpyのarrayとして取り出すことができます。

```Python
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](https://en.wikipedia.org/wiki/Anatomical_plane)

最大値を使うにせよ、最小値を使うにせよ、ここからはnumpy arrayの使いこなしになります。

```Python
# 前略
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画像をなんとかする](http://qiita.com/fukuit/items/ed163f9b566baf3a6c3f)の続きを書いてみました。