病院で学術用途で匿名化された画像ファイルの出力を依頼すると、
/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))
