Help us understand the problem. What is going on with this article?

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

More than 3 years have passed since last update.

DICOM画像とは、簡単にいうと、医療用の画像を装置間でやりとりするための通信プロトコル(DICOM)上で扱われる画像のことで、画像にいろんな付加情報(患者さんの情報であったり、装置の種類であったり、その画像そのものの情報であったり、いろいろなメタデータ)をつけたもののことである。

お医者さんで、CTやMRIなどの検査をした時に、病院によってはお願いすると画像をデータとしてくれることがある。この時のフォーマットが拡張子.dcmだったらDICOM画像だ。

PythonでDICOMを扱うための準備

PythonでDICOM画像を扱いたい場合は、pydicomを使うことになるようだ。

% anaconda show conda-forge/pydicom
Using Anaconda API: https://api.anaconda.org
Name:    pydicom
Summary: Pure python package for DICOM medical file reading and writing
Access:  public
Package Types:  conda
Versions:
   + 0.9.9

To install this package with conda run:
     conda install --channel https://conda.anaconda.org/conda-forge pydicom

というワケなので、conda install -c https://conda.anaconda.org/conda-forge pydicom する。

DICOM画像のサンプルを用意しておく

自分のMRI画像を残念ながら持っていない場合は、sample画像をダウンロードしておく。
DICOM Viewerとして有名なOsirixのサイトに、DICOM Image Libraryがあるので、ここから適当なヤツをもらっておく。
とりあえず、BRAINIXという脳腫瘍のMRI像をダウンロードしてみた。
この.zipファイルを展開すると、MRIで撮像後に、読影医が見るためにいろいろ画処理した画像(T1強調とか、そういうヤツ)があるのだけれど、この中で、「T1-3D-FFE-C\ -\ 801/」というフォルダに入っているファイル100枚を使うことにする。

DICOM画像を読む

以下に、pydicomの機能を使って画像を読み込んでみて、あれこれする例を示す。

DICOM画像を1枚読んでメタデータを参照する

画像を読み込んで、DICOMヘッダを表示させる。

import dicom

d = dicom.read_file('BRAINIX/IM-0001-0075.dcm')
print(d)

ちなみに、特定のメタデータを表示させたい時は、次のようにする。

print(d.Manufacturer)

どうやら、装置メーカーは「Philips Medical Systems」らしい。

print(d[0x0018, 0x0050])

どうやら、Slice厚は1.5mmになっているようだ。

20161030_002.png

このように、DICOMのメタデータにはタグと呼ばれるキー名が決まっていて、キーには番号が振られている。そして、pydicomでは、キー名でも番号でも表示させることができる。

DICOM画像を1枚読んで画像を表示する

DICOM画像の画像部分の実体はnumpyのarray形式で取り出すことができるので、opencvでもmatplotlibでも使えば表示することができる。

Kobito.MQX9VW.png

複数のDICOM画像を読みこんで1枚の連続した断層画像に変換する

ところで、今回利用したBRAINIXの画像は1.5mm厚の連続した100枚の画像になっている。
この100枚の画像を読み込んで一つの3次元配列に取り込むことができる。

import dicom
from matplotlib import pyplot as plt
%matplotlib inline

root_dir = './BRAINIX'
dcms = []
for d, s, fl in os.walk(root_dir):
    for fn in fl:
        if ".dcm" in fn.lower():
            dcms.append(os.path.join(d, fn))
ref_dicom = dicom.read_file(dcms[0])
d_array = np.zeros((ref_dicom.Rows, ref_dicom.Columns, len(dcms)), dtype=ref_dicom.pixel_array.dtype)
for dcm in dcms:
    d = dicom.read_file(dcm)
    d_array[:, :, dcms.index(dcm)] = d.pixel_array

要するに、BRAINIXフォルダ以下の拡張子が.dcmになっているファイルを全て読み込んで、順にd_arrayという三次元配列に内容をコピーしている。

任意のスライスを表示したい場合は、例えば50枚目のスライスならば、pyplt.imshow(d_array[:, :, 49])するだけだ。
これで、3次元配列に読み込んだので、3断面表示も可能になる。

import os
import dicom
from matplotlib import pyplot as plt
%matplotlib inline

# BRAINIXフォルダ以下のdcmファイルを読み込む
root_dir = './BRAINIX'
dcms = []
for d, s, fl in os.walk(root_dir):
    for fn in fl:
        if ".dcm" in fn.lower():
            dcms.append(os.path.join(d, fn))
ref_dicom = dicom.read_file(dcms[0])
d_array = np.zeros((ref_dicom.Rows, ref_dicom.Columns, len(dcms)), dtype=ref_dicom.pixel_array.dtype)
for dcm in dcms:
    d = dicom.read_file(dcm)
    d_array[:, :, dcms.index(dcm)] = d.pixel_array

# 3断面表示する
print(d_array.shape)
plt.subplot(1, 3, 1)
plt.imshow(d_array[127, :, :])
plt.subplot(1, 3, 2)
plt.imshow(d_array[:, 127, :])
plt.subplot(1, 3, 3)
plt.imshow(d_array[:, :, 49])

20161030_003.png

DICOM画像を加工して保存する

読み込んだDICOM画像を加工してもとのDICOM画像に上書き保存する場合は、pixel_array の配列を編集することになる。

ただし、pixel_arrayの配列を編集するだけではもとの画像は変わらない。これは、DICOM画像のイメージデータの実体はPixelDataという配列に保存されており、その読み出し用のインターフェイスがpixel_arrayになっているという事情による。そこで、編集したpixel_arrayの内容が例えばimgというアレイに保存されているのであれば、ソレをPixelDataに書き戻す処理が必要になる。

d = dicom.read_file(dcm)
img = d.pixel_array
# imgの加工
d.PixelData = img.tostring()
d.save_as(dcm)

DICOM画像を円形に切り抜いてみる

「DICOM画像を中央部分のみ、画像の縦横の長さの半分の直径の円で切り抜く」という作業を、何のためか分からないけれど、やってみることにする。

以下のコードで、画像の縦横の半分の直径の円の内側が1で、外側が0の画像ができる。
半径rとすると、画像の縦横サイズはref_dicom.pixel_array.shapeから取得できるので、そのサイズの1/4を設定する。

r = ref_dicom.pixel_array.shape[0] / 4
x, y = np.indices((ref_dicom.Rows, ref_dicom.Columns))
circle = (x - (ref_dicom.Columns / 2))**2 + (y - (ref_dicom.Rows / 2))**2 < r**2
mask = circle.astype(int)

このmaskを掛け算すると、円形の内側が1で、外側が0なので、円形に切り抜くことができる。

### 複数のDICOM画像ファイルを読み込んで、中央部を円形に切り抜く

import os
import dicom
from matplotlib import pyplot as plt
%matplotlib inline

# BRAINIXフォルダ以下のdcmファイルを読み込む
root_dir = './BRAINIX'
dcms = []
for d, s, fl in os.walk(root_dir):
    for fn in fl:
        if ".dcm" in fn.lower():
            dcms.append(os.path.join(d, fn))
ref_dicom = dicom.read_file(dcms[0])

# 円に切り抜くためのマスクを作成する
r = ref_dicom.pixel_array.shape[0] / 4
x, y = np.indices((ref_dicom.Rows, ref_dicom.Columns))
circle = (x - (ref_dicom.Columns / 2))**2 + (y - (ref_dicom.Rows / 2))**2 < r**2
mask = circle.astype(int)

# 全ての画像にmask処理を実施する
d_array = np.zeros((ref_dicom.Rows, ref_dicom.Columns, len(dcms)), dtype=ref_dicom.pixel_array.dtype)
for dcm in dcms:
    d = dicom.read_file(dcm)
    img = d.pixel_array * mask
    d_array[:, :, dcms.index(dcm)] = img

# 3断面表示する
plt.figure(figsize=(8, 6))
plt.subplot(1, 3, 1)
plt.imshow(d_array[127, :, :])
plt.subplot(1, 3, 2)
plt.imshow(d_array[:, 127, :])
plt.subplot(1, 3, 3)
plt.imshow(d_array[:, :, 49])

20161030_004.png


PythonでDICOM画像を扱う時は、だいたいこんな感じ。

本日のコード

fukuit
最近、事務系の職場に異動したので、職業プログラマではなくなりました。
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away