tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート へ

dicomファイルをまとめて変換するスクリプト

dicom opencv 画像処理 matplotlib

病院で学術用途で匿名化された画像ファイルの出力を依頼すると、

/run/media/tak0kada/dvd/DICOM
|--00000000
|  |--00000000
|  |--00000001
|--00000001
...
|--_DCM_INF
|  |--00000000
|  |  |--_DCM_IDX.txt
|  |  |--overlay.txt
...

というディレクトリ構成のDVDをもらえる。dicomファイルは00000000、00000001と連番になっており、画像/動画データである。スクリプトでまとめて変換したい。pydicomだと難しそうだったがsitkとopencvで動画も出力できた。以下詳細。

SimpleITK

病院にあるdicomファイルを吐き出す機械として、CT、MRI、PET-CT、エコーが挙げられる。CTは肺野条件、縦隔条件の2種類のウィンドウを含む白黒の画像、MRIは白黒の画像を出力する。PET-CTはPETだけの白黒画像、CTだけの白黒画像、合成されたカラー画像を出力する。エコーはカラーの動画、白黒あるいはカラーの画像を出力するはずである(エコー室のエコーは高性能なためカラー動画しか出力しない)。この仮定の元、条件分岐はベタ書きした。StateパターンやStrategyパターンで書き直すべきかもしれない。

Pydicomを使ってもデータは取り込めるが、(3,n,m)次元のnumpy.ndarrayが返るため、matploblibに渡す(n,m,3)次元の行列を作り直さないといけないのと、どうも出力がおかしかったことからSimpleITKを使うことにした。Pydicomでも白黒画像は普通に使いまわせたが、import dicomだけでなく、

# 参考: https://groups.google.com/forum/#!topic/pydicom/ffuPvuRQq0s
import dicom.charset
dicom.charset.python_encoding['ISO_IR 13']='cp932' # cp932の代わりにshift-jisでも問題なさそう

としないとファイルの読み込みに失敗した。その後は、ds=dicom.read_file('...')、ds.pixel_arrayとアクセスできる。

SimpleITKを使ってなんとかするのが以下のスクリプト

import cv2
import SimpleITK as sitk
from matplotlib import pyplot as plt
import numpy as np
from pathlib import Path
import os


INPUT_DIR = Path('/run/media/username/dvd/DICOM/')
OUTPUT_DIR = Path('~/dicom')
OUTPUT_DIR.mkdir()


for pp in INPUT_DIR.iterdir():
    if pp.name == 'DCM_INF' or pp.is_file():
        continue
    else:
        Path(OUTPUT_DIR, pp.name).mkdir()

    for p in pp.iterdir():
        os.chdir(str(Path(OUTPUT_DIR, pp.name)))
        data = sitk.ReadImage(str(p))
        img = sitk.GetArrayFromImage(data)


        # モノクロ画像
        if img.shape[0] == 1 and img.ndim == 3:
            # Photometric Interpretationについてはhttps://www.dabsoft.ch/dicom/3/C.7.6.3.1.2/を参照
            cmap = 'binary' if 'MONOCHROME1' in data.GetMetaData('0028|0004') else 'binary_r'

            # ウィンドウ幅についてはhttps://www.escape.gr/manuals/qmedical/commands/index_Values_of_Interest__.html
            # window, centerについて2つの値(縦隔条件、肺野条件)が設定される場合がある:
            #   https://groups.google.com/forum/#!topic/comp.protocols.dicom/xANzkKaD0b
            if '\\' in data.GetMetaData('0028|1050'):
                window_center = [float(i) for i in data.GetMetaData('0028|1050').split('\\')]
                window_width = [float(i) for i in data.GetMetaData('0028|1051').split('\\')]
                # 2種類のウィンドウについてそれぞれ出力
                for i in range(2):
                    # ピクセルの値の範囲を設定しないと正規化されてしまう
                    fig = plt.imshow(
                            img[0], vmin=window_center[i]-window_width[i]/2, vmax=window_center[i]+window_width[i]/2,
                            cmap=cmap, interpolation='nearest'
                            )
                    plt.axis('off')
                    fig.axes.get_xaxis().set_visible(False)
                    fig.axes.get_yaxis().set_visible(False)
                    # 縦隔条件のファイル名にa、肺野条件のファイル名にbを追加
                    plt.savefig('ab'[i] + p.name + '.png', bbox_inches='tight', pad_inches=-0.05)
                    plt.clf()

            # ウィンドウが1種類の場合
            else:
                window_center = float(data.GetMetaData('0028|1050'))
                window_width = float(data.GetMetaData('0028|1051'))
                # ピクセルの値の範囲を設定しないと正規化されてしまう
                fig = plt.imshow(
                        img[0], vmin=window_center-window_width/2, vmax=window_center+window_width/2,
                        cmap=cmap, interpolation='nearest'
                        )
                plt.axis('off')
                fig.axes.get_xaxis().set_visible(False)
                fig.axes.get_yaxis().set_visible(False)
                plt.savefig(p.name + '.png', bbox_inches='tight', pad_inches=-0.05)
                plt.clf()


        # カラー画像は普通のRGBと仮定
        elif img.shape[0] == 1 and img.ndim == 4 and ('RGB' in data.GetMetaData('0028|0004')):
            fig = plt.imshow(img[0], interpolation='nearest')
            plt.axis('off')
            fig.axes.get_xaxis().set_visible(False)
            fig.axes.get_yaxis().set_visible(False)
            plt.savefig(p.name + '.png', bbox_inches='tight', pad_inches=-0.05)
            plt.clf()


        # モノクロ動画
        elif img.shape[0] > 1 and img.ndim == 3:
            raise NotImplementedError()


        # カラー動画
        # 参考: https://github.com/ContinuumIO/anaconda-issues/issues/223
        elif img.shape[0] > 1 and img.ndim == 4 and ('RGB' in data.GetMetaData('0028|0004')):
            frame_rate = int(data.GetMetaData('0008|2144'))
            fourcc = cv2.VideoWriter_fourcc(*'MJPG')
            video = cv2.VideoWriter(p.name + '.avi', fourcc, frame_rate, img.shape[1:3][::-1])

            for i in range(img.shape[0]):
                video.write(img[i])
            video.release()

YAKAMI DICOM ツール集

何とかしようとするもコマンドラインでうまく動かず挫折したスクリプトが以下。左記のリンク先と同様のコマンドだけ独立して実行したところ、エラーを吐かずに終了するものの何も出力されなかったので、pythonの部分にバグがないかどうかすら未確認。

# coding: cp932
# 小さいファイルを集めて動画になっているもの用のスクリプト

import os

# 入出力フォルダ(手抜きのため決め打ち)
OUTPUT_ROOT = u"D:/Users/usr_name/Desktop/output"
INPUT_ROOT = "E:/DICOM"

# ファイルの選択
files = os.listdir(INPUT_ROOT)

print(u"このスクリプトは画像ファイルを集めて動画のように見せるタイプのdicomファイルを変換するものです。")
print(u"ファイルを以下から選択して下さい")

files = [i for i in os.listdir(INPUT_ROOT) if os.path.isdir(INPUT_ROOT + "/" + i)]
for i, dir in files:
    print(i, ": ", dir)
choice = os.listdir([int(input())])

if os.path.exists(OUTPUT_ROOT):
   pass
else:
   os.mkdir(OUTPUT_ROOT)

OUTPUT_DIR = OUTPUT_ROOT + "/" + choice
INPUT_DIR = INPUT_ROOT + "/" + choice
os.mkdir(OUTPUT_DIR)

# DICOMファイルの変換
COMMAND_TEMPLATE = """\
\"C:\\Program Files\\YAKAMI Software\\YAKAMI DICOM Tools (x64 ja)\\cmd_conv.exe\" ^
/mode file_file ^
/infile \"E:\\DICOM\\{choice}\\{file}\" ^
/outfile \"D:\\Users\\usr_name\\Desktop\\output\\{choice}\\{file}\" ^
/fmt png
"""

dicom_pics = [i for i in os.listdir(INPUT_DIR) if not os.path.isdir(INPUT_DIR + "/" + i)]

for file in dicom_pics:
    os.system(COMMAND_TEMPLATE.format(choice = choice, file = file))