2022年10月19日水曜日

Pyゼミ0.06　備忘録　Ubuntu20.04にCUDA11.4、TensorFlowをインストールする　

はじめに

Ubuntuを20.04にして、NVIDIAのドライバーをインストールするとCUDA12.0がインストールされます(2023年2月20日）。

問題なくインストールされnvidia-smiを正常に起動されます。

しかし、TensorflowがGPUを認識しません。

そこで、今回CUDAを11.4にバージョンダウンして利用する方法を投稿します。

1）Ubuntu20.04インストール後　ソフトウェアアップデート

$ sudo apt update
$ sudo apt upgrade

$ sudo reboot

2)CUDA＆Driverインストール

下記URLにてCUDA Toolkit Archiveのページへ
https://developer.nvidia.com/cuda-toolkit-archive

バージョン11.4の下記を選択する

CUDA Toolkit Documentation v11.4.4

Select Target Platformは以下を選択

OS：Linux
Architecture：x86_64
Distribution：Ubuntu
Version：20.04
Installer Type：dev(local)

注意：cudaのインストール時に明示的に11.4をcuda-11-4と指定すること

$ wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/cuda-ubuntu2004.pin
$ sudo mv cuda-ubuntu2004.pin /etc/apt/preferences.d/cuda-repository-pin-600
$ wget https://developer.download.nvidia.com/compute/cuda/11.4.4/local_installers/cuda-repo-ubuntu2004-11-4-local_11.4.4-470.82.01-1_amd64.deb
$ sudo dpkg -i cuda-repo-ubuntu2004-11-4-local_11.4.4-470.82.01-1_amd64.deb
$ sudo apt-key add /var/cuda-repo-ubuntu2004-11-4-local/7fa2af80.pub
$ sudo apt-get update
$ sudo apt-get -y install cuda-11-4

$ sudo reboot

確認１

$ nvidia-smi

確認２

/usr/localにcudaの関連ディレクトリがあるか確認する

$ sudo find /usr/local | grep cuda

/usr/local/cuda-11.4
/usr/local/cuda-11.4/compute-sanitizer
・・・

3)CUDAパス設定

CUDA のライブラリにパスを通します

$ gedit .bashrc

次の３行を　.bashrc　ファイルの最後に追加

#CUDA
export PATH="/usr/local/cuda/bin:$PATH"
export LD_LIBRARY_PATH="/usr/local/cuda/lib64:$LD_LIBRARY_PATH"

確認

$ source .bashrc
$ nvcc -V

nvcc: NVIDIA (R) Cuda compiler driver
Copyright (c) 2005-2021 NVIDIA Corporation
Built on Mon_Oct_11_21:27:02_PDT_2021
Cuda compilation tools, release 11.4, V11.4.152
Build cuda_11.4.r11.4/compiler.30521435_0

4)cuDNNのインストール

下記のURLからNVIDIA Developer Program Membership Requiredへ
https://developer.nvidia.com/rdp/cudnn-download

ログイン（またはJoin nowにて新規登録）

cuDNN Downloadのダウンロードページへ

今回はcuDNN v8.6をインストールするため、最新のボタンの下の下記を選択
Archived cuDNN Releases

リストから下記を選択しファイルをダウンロード

Download cuDNN v8.6.0 (October 3rd, 2022), for CUDA 11.x
Local Installer for Ubuntu20.04 x86_64 (Deb)

ダウンロードしたファイル名：cudnn-local-repo-ubuntu2004-8.6.0.163_1.0-1_amd64.deb

ファイルをhomeディレクトリに移動し、以下を実行

$ sudo dpkg -i cudnn-local-repo-ubuntu2004-8.6.0.163_1.0-1_amd64.deb

$ sudo cp /var/cudnn-local-repo-ubuntu2004-8.6.0.163/cudnn-local-B0FE0A41-keyring.gpg /usr/share/keyrings/
$ sudo apt update
$ sudo apt -y install libcudnn8 libcudnn8-dev libcudnn8-samples

確認

$ dpkg -l | grep cudnn

ii cudnn-local-repo-ubuntu2004-8.6.0.163 1.0-1 amd64 cudnn-local repository configuration files
ii libcudnn8 8.6.0.163-1+cuda11.8 amd64 cuDNN runtime libraries
ii libcudnn8-dev 8.6.0.163-1+cuda11.8 amd64 cuDNN development libraries and headers
ii libcudnn8-samples 8.6.0.163-1+cuda11.8 amd64 cuDNN samples

5）Anacondaのインストール

Python 3.9をインストールするため、昨年のファイルを使用

https://www.anaconda.com/distribution/

ダウントード：Anaconda3-2023.03-Linux-x86_64.sh

$ bash Anaconda3-2023.03-Linux-x86_64.sh

確認

$ source .bashrc

$ conda -V
conda 4.12.0

$ python
[GCC 7.5.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>>

Python 3.9.12 (main, Apr 5 2022, 06:56:58)

6）TnesorFlow と Kerasのインストール

Tensorflowの2.8.2をインストール。

次を実行することでGPU対応のTensorflowとKerasがインストールされる

$ pip install tensorflow==2.8.2

GPU確認

python起動後

>>> from tensorflow.python.client import device_lib
>>> device_lib.list_local_devices()

2023-03-26 14:06:15.580531: I tensorflow/core/platform/cpu_feature_guard.cc:151] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations: AVX2 FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.
2023-03-26 14:06:15.612622: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:15.706883: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:15.708182: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:16.083901: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:16.084200: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:16.084453: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:16.084689: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1525] Created device /device:GPU:0 with 5239 MB memory: -> device: 0, name: NVIDIA GeForce GTX 1060 6GB, pci bus id: 0000:01:00.0, compute capability: 6.1
[name: "/device:CPU:0"
device_type: "CPU"
memory_limit: 268435456
locality {
}
incarnation: 3029067233870127766
xla_global_id: -1
, name: "/device:GPU:0"
device_type: "GPU"
memory_limit: 5493751808
locality {
bus_id: 1
links {
}
}
incarnation: 7369121378162548002
physical_device_desc: "device: 0, name: NVIDIA GeForce GTX 1060 6GB, pci bus id: 0000:01:00.0, compute capability: 6.1"
xla_global_id: 416903419
]
>>>

7)OpenCVのインストール

$ pip install opencv-python

Collecting opencv-python
Downloading opencv_python-4.6.0.66-cp36-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (60.9 MB)
|████████████████████████████████| 60.9 MB 10.5 MB/s
Requirement already satisfied: numpy>=1.14.5 in ./anaconda3/lib/python3.9/site-packages (from opencv-python) (1.21.5)
Installing collected packages: opencv-python
Successfully installed opencv-python-4.6.0.66

確認

>>> import cv2

>>> cv2.__version__
'4.7.0'

8)Pydicomインストール

$ pip install pydicom

Collecting pydicom
Downloading pydicom-2.3.0-py3-none-any.whl (2.0 MB)
|████████████████████████████████| 2.0 MB 3.1 MB/s
Installing collected packages: pydicom
Successfully installed pydicom-2.3.1

9)DCMTKのインストール

$ sudo apt install dcmtk

確認

$ dcmdjpeg --version

（以下の内容が表示されればOK）

$dcmtk: dcmdjpeg v3.6.4 2018-11-29 $
dcmdjpeg: Decode JPEG-compressed DICOM file
Host type: Debian
Character encoding: UTF-8
External libraries used:
- ZLIB, Version 1.2.11
- IJG, Version 6b 27-Mar-1998 (modified)

10)ImageJのインストール

ImageJのページからLinux版をダウンロード

https://imagej.nih.gov/ij/

ダウンロード：ij153-linux64-java8.zip

Homeで展開、起動してアイコンをお気に入りに登録する

GUI画面が小さい場合

Edit/Options/AppearanceのGUI scaleを1.5にしてみる。

11）その他の設定

○メニューに空のドキュメントを追加する

# フォルダーが日本語のままの場合

$ touch ~/テンプレート/空のドキュメント

○エディタ

・行番号の表示

・Tab　4文字空白で設定

・フォントサイズ　16

-----------------------------------

おまけ）Dropboxのインストール

$ cd ~ && wget -O - "https://www.dropbox.com/download?plat=lnx.x86_64" | tar xzf -

起動

$ ~/.dropbox-dist/dropboxd

2022年5月6日金曜日

Pyゼミ1.08　PET画像を読む

医療画像の画素値を定量評価に使う

PET画像など画素値を定量評価に使う場合があります。

以前に紹介した「Pyゼミ1.02　DICOMからJPEGなどへの変換」では、CT画像をウインドニングして0から255の値に変換しています。

したがって、本来のCT値を読むことはできません。

今回はPET画像を対象にPETの本来の画素値を取り出す方法を紹介します。

※PET, Positron Emission Tomography（陽電子放出断層撮影）

プログラミング　スタート

プログラムはここから

PET画像はここから

初めに準備です。

importするライブラリは以下の通りです。

本記事ではpydicomのapply_modality_lutの機能を紹介しますのでimportしておきます。

ファイル名を変数dcmfnmにセットしておきます。

import pydicom

import numpy as np

import cv2

from pydicom.pixel_data_handlers.util import apply_modality_lut

dcmfnm = "./an_pet.50.dcm"

PET画像を読み16ビットの画素値を取り出す

PETの1画素に割り当てらえているビット数は16ビット（Bits Allocated)です。

そして有効なデータとして15ビット（Bits Stored)を使って格納されています。

16ビットの深さでリスケール変換して取り出すメソッドreadimageWindowing16メソッドを作成します。

def readimageWindowing16(filepath):
print("Image Widowing16")
ds = pydicom.dcmread(filepath)
#print("ds=\n",ds)
img = ds.pixel_array
print("\timg shape=",img.shape, type(img))
print("\tOrg img max=",np.max(img),"min=",np.min(img))

ri = ds.RescaleIntercept
rs = ds.RescaleSlope

img = img*rs + ri # Convert pixel value to PET value
img = img.astype(np.uint16) # Unsigned 16bit
print("\trsIntersept =",ri, "\n\trsSlope =",rs)

return img

img16 = readimageWindowing16(dcmfnm)

print("\nWindowing 16bit img: max =", np.max(img16), "min=",np.min(img16), "type=",img16.dtype)

cv2.imwrite("win16.png", img16)

このメソッドは、引数filepathにDICOMファイルのpathを与えると、PETの正しい画素値に変換された画像データが返ります。

dcmreadメソッドを使い、DICOMのデータセットをdsに読み込んでいます。

そして画像データpixel_arrayをimgに代入しています。

正しいPETの画素値に変換するためにリスケール切片とリスケール傾斜をそれぞれ変数riとrsに代入します。

PETの画素値への変換は、読み込んだ画像imgにリスケール傾斜rsを掛けてリスケール切片riを加えるだけです。

最後にこのメソッドを呼び出し、PETの正しい画素値に変換された画像をimg16に代入しています。

また、変換後の最大と最小を表示して、PNGファイル（win16.png)に保存しています。

実行してみましょう

実行するとコンソールに以下のように表示されます。

Image Widowing16

img shape= (128, 128) <class 'numpy.ndarray'>

Org img max= 32767 min= 0

rsIntersept = 0

rsSlope = 1.728227

Converted Windowing 16bit img: max = 56628 min= 0 type= uint16

画像サイズは128×128です。

最大画素値は32,767です。これは2進数で見ると15ビットすべて1になります。

つまり15ビット長の最大値です。

最小画素値は0です。

変換後の最大画素値（PETの画素値）は56,628となっています。

これは画素値にリスケール傾斜を掛けているために値が1.728227倍になっています。

画素値を格納するビット数が15ビットの場合、最大値は32,767に決定されます。

しかし、実際には15ビットで扱える最大値を超える場合があります。このとき、リスケール傾斜で調整を行っています。

結果画像を見てみましょう。

PETファントムで収集した再構成画像です。きれいに表示されています。

img16の配列の各要素にはPETの各値が格納されています。

この値を使って定量評価ができますね。

apply_modality_lutを使ってみる

pydicomのpixel_data_handlers.utilライブラリの中にapply_modality_lutメソッドがあります。

これは、読み込んだ画像に対して簡単にモダリティのLUTやリスケール操作を行うことができます。

それではreadimageModalityLutメソッドを作成して、このメソッドを呼び出してます。

変換された画像imgModの最大値と最小値を表示して、画像を保存します。

def readimageModalityLut(filepath):

print("\nConvert with Modalyty_Lut")

ds = pydicom.dcmread(filepath)

pxarry = ds.pixel_array

img = apply_modality_lut(pxarry, ds)

img = img.astype(np.uint16)

return img

imgMod = readimageModalityLut(dcmfnm)

print("\nConverted ModalityLut img: max =", np.max(imgMod), "min=",np.min(imgMod), "type=",imgMod.dtype)

cv2.imwrite("petModLut.png", imgMod)

実行してみましょう

実行すると次のようにコンソール表示されます。

Convert with Modalyty_Lut

Converted ModalityLut img: max = 56628 min= 0 type= uint16

変換後の最大値は56,628となっています。

先のreadimageWindowing16メソッドの結果と同じです。

しかし、プログラムはリスケール傾斜や切片を読むことなく画素データpxarryをapply_modality_lutメソッドで直接PETの画素値に変換しています。

画像も以下のようにきれいに表示されています。

8ビットでウインドニングしてしまう

16ビットの画素値を取り出さなくてはいけないデータを8ビットでウインドニングして出力するとどうなるでしょう。

まず、8ビットでウインドニングするメソッドを作成します。

def readimageWindowing8(filepath):

print("Image Widowing 8bit")

ds = pydicom.dcmread(filepath)

#print("ds=\n",ds)

img = ds.pixel_array

print("\timg shape=",img.shape, type(img))

print("\tOrg img max=",np.max(img),"min=",np.min(img))

wc = ds.WindowCenter

ww = ds.WindowWidth

ri = ds.RescaleIntercept

rs = ds.RescaleSlope

img = img*rs + ri # Convert pixel value to PET value

print("\twc=",wc,type(wc), "\n\tww=", ww, type(ww))

print("\trsIntersept =",ri, "\n\trsSlope =",rs)

maxval = wc + ww//2

minval = wc - ww//2

img8 = (img - minval)/(maxval - minval)*255

img8 = np.clip(img8, 0, 255).astype(np.uint8)

return img8

img8 = readimageWindowing8(dcmfnm)

print("\nConverted Windowing 8bit img: max =", np.max(img8), "min=",np.min(img8), "type=",img8.dtype)

cv2.imwrite("win8.png", img8)

readimageWindowing8メソッドの中で画像imgをリスケール傾斜と切片を使って変化するところまでは同じです。

その後、ウインドウセンタwcとウインドウ幅wwを用いて表示する最大値maxvalと最小値minvalを決めています。

そして画像のウインドニング処理を行っています。

実行してみると次のようにコンソールに表示されます。

Image Widowing 8bit

img shape= (128, 128) <class 'numpy.ndarray'>

Org img max= 32767 min= 0

wc= 6351.233 <class 'pydicom.valuerep.DSfloat'>

ww= 12702.47 <class 'pydicom.valuerep.DSfloat'>

rsIntersept = 0

rsSlope = 1.728227

Converted Windowing 8bit img: max = 255 min= 0 type= uint8

ウインドウセンタwcとウインドウ幅がそれぞれ6351.233と12702.47となり、これから計算される最大値maxvalと最小値minvalは12702.233と0.2330となります。最大値は変換後の画素の持つ最大値56,628に遠く及びません。

出力したPNG画像は白くつぶれた画像になってしまいます。

おわりに

ある先生からPETの定量解析をしたいという要望があり本記事を紹介しました。

以前紹介した8ビット画像で出力する方法ですと最大値は255に固定されてしまいます。

また、ウインドセンタとウインドウ幅の適切な値がDICOMタグに入っていない場合出力した画像自体おかしな表示になります。

ImageJでも表示しましたがもちろん同じ結果でした。

この原因として考えられるのは、最初に表示した画像でウインドウセンタとウインド幅を調整した後、それ以降の画像は調整されず、過去の調整結果が引き継がれたのでしょう。

ファントム実験のデータなのでそういうこともありそう。

最大値が56,628、最小値が0の画像を表示するのであれば、適切なウインドウセンタが28,314でウインドウ幅が56,628となるはずです。

ブログに上げようと思い、掲載まで約1年かかってしまいました。

ほんとうに年月の経つのは早いですね。

2022年5月2日月曜日

Pyゼミ1.0　番外編　DICOM RDSRを読む

Pydicomを使ってRDSRを読んで必要な情報を取り出してみる

RDSRはRadiation Dose Structured Reportの略で、CTや血管造影などX線を用いた検査でどのくらX線量（Radiation dose）を照射したのか、構造化された書式（Structured report）にて出力したものです。

RDSRはDICOM規格で規定されています。詳細は規格書を読みましょう。

今回は時間がなく規格書を理解せず、帰納的にSRの解析に取り組みました。

解析結果が正しいかどうか、医療現場で使用しているソフトウェアの結果で検証しています。

検証した範囲では問題なさそうですが、利用は自己責任でお願いします。

プログラムはこちらから

RDSRファイルはこちらから

（この記事は2022年4月24日に日本放射線技術学会北海道支部第78回春季学術大会　医療情報セミナーで発表したものです）

PydicomはDICOM画像だけでなくSRも読むことができます。RDSRファイルを読み込んで表示するだけであれば数行で書くことができます。

import pydicom
fnm = "./CT_rdsr.SR"
ds = pydicom.dcmread(fnm)
print(ds)

このたった4行でRDSRのDICOMファイル"rdsr.sr"の内容をpydicomを用いて変数dsに読み込み、それをプリントすることでRDSRの内容を確認することができる。

DICOM RDSRの内容を見てみる

DICOMのファイルはデータセット呼ばれるデータエレメントの集合でできています。
データエレメントは、タグと呼ばれる16進数のグループ番号とエレメント番号から始まり、キーワード、値表現（VR）、値の長さ（VL）と値（VF, Value field）で構成されています。

print(ds)で出力された内容を見てみしょう。

表示された各行はデータエレメントを表していて次のように表示されています。

(グループ番号, エレメント番号) キーワード VR:値

(0002,0002)と(0008,0016)のSOP Class UIDが”X-Ray Radiation Dose SR Storage"と記述されています。

また、(0008,0060)のModalityは”SR”とあり、このファイルがStructured reportであることがわかります。

Dataset.file_meta -------------------------------

(0002, 0000) File Meta Information Group Length UL: 196

(0002, 0001) File Meta Information Version OB: b'\x00\x01'

(0002, 0002) Media Storage SOP Class UID UI: X-Ray Radiation Dose SR Storage

(0002, 0003) Media Storage SOP Instance UID UI: 1.3.12.2.1107.5.1.4.65939.30000022022322414137700000249

(0002, 0010) Transfer Syntax UID UI: Explicit VR Little Endian

(0002, 0012) Implementation Class UID UI: 1.3.12.2.1107.5.4.7

(0002, 0013) Implementation Version Name SH: 'SIEMENSXLEOVB21C'

-------------------------------------------------

(0008, 0005) Specific Character Set CS: 'ISO_IR 100'

(0008, 0016) SOP Class UID UI: X-Ray Radiation Dose SR Storage

(0008, 0018) SOP Instance UID UI: 1.3.12.2.1107.5.1.4.65939.30000022022322414137700000249

(0008, 0020) Study Date DA: '20220224'

(0008, 0021) Series Date DA: '20220224'

(0008, 0023) Content Date DA: '20220224'

(0008, 0030) Study Time TM: '074907.240000'

(0008, 0031) Series Time TM: '075458.275000'

(0008, 0033) Content Time TM: '075458.315000'

(0008, 0050) Accession Number SH: '6010016637053'

(0008, 0060) Modality CS: 'SR'

上記以外の情報は通常のDICOM画像とそう変わりありません。

しかし、もう少し下のほうを見てみると、(0040,0275)のRequest Attributes Sequence以下で字下げ（インデント）が見られ、

さらに(0040,0008) Scheduled Protocol Code Sequence以下でインデントが発生しています。

これはVRが"SQ"に設定されているためです。

”SQ”は「ある意味のある集合のシーケンス（連続した項目）を含み、このシーケンスは入れ子状態をとることができる」ものです。

詳しくはDICOM Part５（http://www.jira-net.or.jp/dicom/dicom_data_02_01.html）を参照してください。

(0029, 1008) [syngo Report Type]                 CS: 'CTDOSEREPORT'

(0029, 1009) [syngo Report]                      LO: '1.0'

(0029, 1015) [SR Variant]                        US: 6

(0029, 1017) [SC SOP Instance UID]               UI: 1.3.12.2.1107.5.1.4.65939.30000022022322414137700000247

(0040, 0275)  Request Attributes Sequence  1 item(s) ---- 

   (0040, 0008)  Scheduled Protocol Code Sequence  1 item(s) ---- 

      (0008, 0100) Code Value                          SH: '3003100550020000'

      (0008, 0102) Coding Scheme Designator            SH: '99'

      (0008, 0103) Coding Scheme Version               SH: '6010016637053'

      (0008, 0104) Code Meaning                        LO: 'AG?F?~???ETAE?F?S?g'

      ---------

   (0040, 0009) Scheduled Procedure Step ID         SH: '6010016637053'

   (0040, 1001) Requested Procedure ID              SH: '6010016637053'

   ---------

この"SQ"をさらに追いかけてみます。

すると、(0040, a730)　Content Sequence が20項目あると表示されています。

つまり、これ以下に何らかの意味ある項目が20項目あることを示しています。

(0040, a050) Continuity Of Content CS: 'SEPARATE'

(0040, a372) Performed Procedure Code Sequence 0 item(s) ----

(0040, a491) Completion Flag CS: 'COMPLETE'

(0040, a493) Verification Flag CS: 'UNVERIFIED'

(0040, a504) Content Template Sequence 1 item(s) ----

(0008, 0105) Mapping Resource CS: 'DCMR'

(0040, db00) Template Identifier CS: '10011'

---------

(0040, a730) Content Sequence 20 item(s) ----

(0040, a010) Relationship Type CS: 'HAS CONCEPT MOD'

(0040, a040) Value Type CS: 'CODE'

(0040, a043) Concept Name Code Sequence 1 item(s) ----

(0008, 0100) Code Value SH: '121058'

(0008, 0102) Coding Scheme Designator SH: 'DCM'

(0008, 0104) Code Meaning LO: 'Procedure reported'

---------

(0040, a168) Concept Code Sequence 1 item(s) ----

(0008, 0100) Code Value SH: 'P5-08000'

(0008, 0102) Coding Scheme Designator SH: 'SRT'

(0008, 0104) Code Meaning LO: 'Computed Tomography X-Ray'

---------

(0040, a730) Content Sequence 1 item(s) ----

(0040, a010) Relationship Type CS: 'HAS CONCEPT MOD'

(0040, a040) Value Type CS: 'CODE'

(0040, a043) Concept Name Code Sequence 1 item(s) ----

(0008, 0100) Code Value SH: 'G-C0E8'

(0008, 0102) Coding Scheme Designator SH: 'SRT'

(0008, 0104) Code Meaning LO: 'Has Intent'

---------

さらに、下を見てみると次のように表示されています。

(0008, 0104) Code Meaningには'CT Accumulated Dose Data'とあり、

その下に

(0008, 0104) Code Meaningには'Total Number of Irradiation Events'とあり、

またその下に

(0008, 0104) Code Meaningには'events'とあり、

次に

(0040, a30a) Numeric Valueには"7.0"とあります。

さらに、その下には

(0008, 0104) Code Meaning には'CT Dose Length Product Total'とあり、

その下に

(0008, 0104) Code Meaningには'mGy.cm'と単位があり、

次に

(0040, a30a) Numeric Valueには"377.94"と数値があります。

(0040, a010) Relationship Type CS: 'CONTAINS'

(0040, a040) Value Type CS: 'CONTAINER'

(0040, a043) Concept Name Code Sequence 1 item(s) ----

(0008, 0100) Code Value SH: '113811'

(0008, 0102) Coding Scheme Designator SH: 'DCM'

(0008, 0104) Code Meaning LO: 'CT Accumulated Dose Data'

---------

(0040, a050) Continuity Of Content CS: 'SEPARATE'

(0040, a730) Content Sequence 2 item(s) ----

(0040, a010) Relationship Type CS: 'CONTAINS'

(0040, a040) Value Type CS: 'NUM'

(0040, a043) Concept Name Code Sequence 1 item(s) ----

(0008, 0100) Code Value SH: '113812'

(0008, 0102) Coding Scheme Designator SH: 'DCM'

(0008, 0104) Code Meaning LO: 'Total Number of Irradiation Events'

---------

(0040, a300) Measured Value Sequence 1 item(s) ----

(0040, 08ea) Measurement Units Code Sequence 1 item(s) ----

(0008, 0100) Code Value SH: '{events}'

(0008, 0102) Coding Scheme Designator SH: 'UCUM'

(0008, 0103) Coding Scheme Version SH: '1.4'

(0008, 0104) Code Meaning LO: 'events'

---------

(0040, a30a) Numeric Value DS: "7.0"

---------

(0040, a010) Relationship Type CS: 'CONTAINS'

(0040, a040) Value Type CS: 'NUM'

(0040, a043) Concept Name Code Sequence 1 item(s) ----

(0008, 0100) Code Value SH: '113813'

(0008, 0102) Coding Scheme Designator SH: 'DCM'

(0008, 0104) Code Meaning LO: 'CT Dose Length Product Total'

---------

(0040, a300) Measured Value Sequence 1 item(s) ----

(0040, 08ea) Measurement Units Code Sequence 1 item(s) ----

(0008, 0100) Code Value SH: 'mGy.cm'

(0008, 0102) Coding Scheme Designator SH: 'UCUM'

(0008, 0103) Coding Scheme Version SH: '1.4'

(0008, 0104) Code Meaning LO: 'mGy.cm'

---------

(0040, a30a) Numeric Value DS: "377.94"

---------

このように、Code Meaningの値やNumeric Valueの値を抽出することで必要な情報を得ることができるのではないかと考えました。

少し解析のフローをまとめると

データエレメントのKeywordが“Content Sequence”が見つかったら解析を始める
データエレメントのVRが“SQ“のとき、KeywordとValueを取得する
データエレメントのkeywordが“Code Meaning”のとき、Valueの値を出力する
取得したValue（データエレメント）について２.の解析を繰り返す。

つまり、VRに"SQ"が見つかったらその下の階の項目を取り出し、その中のVRに"SQ"が見つかったら、さらにその下の。。。と繰り返すことで各階のCode Meaningの情報を取り出すこと考えました。

しかし、階層が3階層より下が上手く取り出すことができませんでした。

これはPydicomの仕様（バグ？）かわかりません。

そこで、タグ、キーワード、VRと値を取り出すメソッド（getInfoDCM）を作成して対応しました。

RDSRを読むPythonプログラム

RDSRを読み込むための準備を最初に書きます。

インポートするのPydicomだけです。

RDSRのファイル名を変数fnmに代入します。

rstxt, dstxtとsmryはSRの処理結果を格納する変数です。

nLvは"SQ"による階層の深さ、インデントの深さを代入する変数です。

prevLは前の"SQ"の行（データエレメント）の深さを格納する変数です。

inkeywordsには、抽出すべき値のキーワードを格納しています。

import pydicom

fnm = "./CT_rdsr.SR" # RDSR DICOM file name

rstxt = "" # For result text using getInfoDCM method

dstxt = "" # For DICOM dataset

smry = "" # For target text

nLv = 0 # Depth Level

prevL = 0 # For depth level of previous SQ line

# For extraction keyword

inkeywords = "CodeMeaning,TextValue,NumericValue,UID,DateTime"

データエレメントの深さ求めるメソッドとデータエレメントの各要素求めるメソッドについて

１．データエレメントの深さを求めるgetDepthメソッド

　インデントは各行（データエレメント）の前の空白の数を数えて返します。

　引数のdataはひとつのデータエレメントです。

def getDepth(data):

cnt = 0

for i in range(len(data)):

if data[i] == ' ':

cnt += 1

else:

break

return cnt

２．データエレメントの各要素求めるgetInfoDCMメソッド

　このメソッドはデータエレメント中のタグ、キーワード、VRと値を抽出して返します。

　引数のdataには　(grp,elm) Key word VR: 'value' のようにひとつのデータエレメントが格納されています。

Pydicomでは"SQ"のシーケンスの終わりに'---------'が表示されます。

メソッドの最初にこの処理を記述します。

与えれたdataの内容が'---------'であった場合、

tagとVRには何も値を入れず、

キーワードには、”NoKeyword”と

値には、'---------'をそのまま入れて返します。

def getInfoDCM(data):

# ex) data is "(grp,elm) Key word VR: 'value' "

data = data.strip(" ") # Remove blanks on both sides of string

if data == '---------':

tag = ["",""]

kwd = "NoKeyword"

vr = ""

val = '---------'

return tag, kwd, vr, val

# Tag

tag = data.split(")")[0] # "(grp, elm"

tag = tag.replace("(", "").replace(" ","") # "grp, elm"

tag = tag.split(",") # [group, elemnt]

# vr

vr = ""

if "item(s)" in data: # ex) (0040, a043) Concept Name Code Sequence 1 item(s) ----

vr = "SQ"

else:

vr = data.split(":")[0][-2:] # ex) "(grp,elm) Key word VR"->"VR"

# value

val = ""

if ":" in data:

val = data.split(":")[-1] # [・・・ VR:value] -> value

elif "item(s)" in data:

val = ' '.join(data.split(" ")[-3:-1]) # -> n item(s)

val = val.replace("\'","").replace("\"","") # remove " and '

# keyword

kwd = ""

stat = data.find(")")+1 # Start position in keyword

end = data.find(vr) # End position in keyword

kwd = data[stat:end] # Extract keyword (remove (g,e))

if "item(s)" in kwd: # Case of VR="SQ"

kwd = kwd.split("item(s)")[0] # [kwd n,----] -> (g,e) kwd n

kwd = kwd.strip(" ").split(" ")[0:-1] # Extract the string before "n items"

kwd = ''.join(kwd).strip(" ") # Joint keyword and remove blanks

else:

kwd = kwd.replace(" ","") # Remove blank in keyword

return tag, kwd, vr, val

はじめにTagの処理を見ます。

(grp,elm) Key word VR: 'value'

グループ番号とエレメント番号は()の中にありますのでこれを頼りに２つの番号を抽出し、リストtagに代入します。

次に、VRです。

(grp,elm) Key word VR: 'value'

(grp,elm) Key word 99 item(s) ----

変数dataの中に"item(s)"があればVRは”SQ”です。

"SQ"以外の場合（else)、コロン（：）で変数dataの値を２つ分割し、前半[0]の最後の2文字[-2:]を抽出し変数vrに代入します。

次に、値です。

(grp,elm) Key word VR: 'value'

(grp,elm) Key word 99 item(s) ----

変数data中にコロン（：）があれば2つに分割し後半[1]を値として変数valに代入します。

コロンがなく、”item(s)”があれば、変数dataを空白で分割し”99 items(s)"のみ抽出しvalに代入します。

最後にvalの値の前後のシングル/ダブルクォーテーション（ ’ / " ）を削除します。

最後にキーワードです。

(grp,elm) Key word VR: 'value'

(grp,elm) Key word 99 item(s) ----

データエレメントの”）”の後ろの文字（stat番目）からVRの前（end)の文字列を抽出data[stat:end]して変数kwdに代入します。

しかし、この処理では、VRが"SQ"のときキーワードに"n item(s)"が存在します。

これを除く処理を行いました。

また、抽出したキーワード中の空白は除くようにしました。つまり”Study Date”は”StudyDate”になります。

これはPydicomのキーワードの表現に揃えるためです。

メインのプログラムを読む

はじめにpydicomを使ってRDSRのファイルを読み込みデータセットを変数dsに代入します。

SOPClassUIDなど必要な基本情報を表示します。

smryは必要な抽出項目の要約を代入する変数です。

変数dsは'pydicom.dataset.FileDataset'の型を取りますが、これを文字列型に変換して変数txdsに代入します。

txdsを改行（\n)で分割して行単位に分割してリストlndsに代入します。

読み込んだ情報などを表示していったん処理が停止（pause）します。

Enterキーで処理が進みます。

for文でlndsから1行（1データエレメント）ずつ取り出して処理します。

valueに'ContentSequence'が見つかったら、smryflagをTrueにして必要情報を抽出するように処理します。

変数rstxtには行番号、深さのレベル数とともに抽出した結果を追加代入(+=)しています。

一方、dstxtにはデータセットに行番号とデータエレメントを追加代入しています。

この２つの出力ファイルを比較することにより、間違いなく抽出できているか確認できます。

変数smryには注目するキーワードのデータエレメントの値を追加代入しています。

# ===== MAIN =====

ds = pydicom.dcmread(fnm) # Read DICOM SR

#print("ds=\n",ds)

print("SOP Class UID:", ds.SOPClassUID)

print("Modality :", ds.Modality)

print("Manufacturer :", ds.Manufacturer)

print("Study Date :", ds.StudyDate)

print("Pt Name :", ds.PatientName)

print("Pt ID :", ds.PatientID)

print("Modality :", ds.Modality)

print("Content Sequence:", len(ds.ContentSequence))

input("pause. Type ENTER key")

smry = "Study Date," + str(ds.StudyDate) + "\n"

smry += "Study Time," + str(ds.StudyTime) + "\n"

txds = str(ds) # Convert pydicom.dataset to String

print(txds)

print("\nData type :",type(ds),"\n\t\t-->",type(txds))

lnds = txds.split("\n") # Split text to lines

print("len(ds)=",len(lnds))

input("pause. Type ENTER key")

for i, ln in enumerate(lnds):

tag, keyword, vr, value = getInfoDCM(ln)

nLv = getDepth(ln)//3

if value in 'ContentSequence': # Start anlysis from this line

smryflag = True

sIndent = " "*(nLv*3) # Indent

rstxt += str(i)+",L:"+str(nLv)+sIndent+"("+str(tag)+"),"+keyword+","+vr+","+value+"\n"

dstxt += str(i)+","+ln + "\n"

print("value = >",value,"<")

# CR flag

if prevL > nLv or nLv <= 1:

crflag = True

else:

crflag = False

if smryflag and keyword in inkeywords:

#smry += "["+str(i)+"L"+str(+nLv)+"]"+str(value)+", "

smry += str(value)+", "

if vr == "SQ" and crflag: # VR is SQ, crflag is True

smry += "\n" # Then add CR in smry

if vr == "SQ": # VR is SQ,

prevL = nLv # Then assign now level to prevL

smry += rmCRline(smry)

with open("rdsrResult.txt","w") as f:

f.write(rstxt)

with open("rdsrDataset.txt","w") as f:

f.write(dstxt)

with open("rdsrSummary.csv","w") as f:

f.write(smry)

さて、ここで問題です。変数smryに値を追加していくだけですと、抽出したい情報がたった1行で出力されてしまいます。

そこで、意味ある情報の塊の境界で改行を入れる必要があります。

改行するかどうかを変数crflagを用いて制御しています。

if prevL > nLv or nLv <= 1:

crflag = True

else:

crflag = False

これは、1つ前の"SQ"のデータエレメントの階層レベルprevLが現在の階層レベルnLvより大きい（つまり現在の階層は前の階層より上がるとき）、または現在の階層が1以下のとき（つまり階層が0または1の深さにある）ときはcrflagをTrueにしています。

それ以外はcrflagをFalseに設定します。

そして、VRの値とcrflagの状態を見て変数smryに改行を入れています。

if vr == "SQ" and crflag: # VR is SQ, crflag is True

smry += "\n" # Then add CR in smry

このままでは改行が至る所に入ってしまうため、無駄な改行を削除するためrmCRlineメソッドを作成しました。

smry += rmCRline(smry)

実行してみる

プログラム実行してみましょう

初めに基本的な情報を表示して停止します。

SOP Class UID: 1.2.840.10008.5.1.4.1.1.88.67

Modality : SR

Manufacturer : SIEMENS

Study Date : 20220224

Pt Name : NO^Name

Pt ID : 12345678

Modality : SR

Content Sequence: 20

pause. Type ENTER key

ここで、「ENTER」キーを打つと、処理が次に進みます。

・・・・・・・

(0008, 0104) Code Meaning LO: 'Source of Dose Information'

---------

(0040, a168) Concept Code Sequence 1 item(s) ----

(0008, 0100) Code Value SH: '113856'

(0008, 0102) Coding Scheme Designator SH: 'DCM'

(0008, 0104) Code Meaning LO: 'Automated Data Collection'

---------

Data type : <class 'pydicom.dataset.FileDataset'>

--> <class 'str'>

len(ds)= 3159

pause. Type ENTER key

読み込んだデータセットのデータの型が'pydicom.dataset.FileDataset'から'str'に変換されました。

「ENTER」キーを打つと、処理が次に進みます。

・・・・・・・

value = > CONTAINS <

value = > CODE <

value = > 1 item(s) <

value = > 113854 <

value = > DCM <

value = > Source of Dose Information <

value = > --------- <

value = > 1 item(s) <

value = > 113856 <

value = > DCM <

value = > Automated Data Collection <

value = > --------- <

最後に３つのファイルが出力されます。

rdsrDataset.txt データエレメントに行番後を付加して出力します。

rdsrResult.txt getInfoDCMメソッドで抽出した結果を行番号と一緒に出力します。

rdsrSummary.txt 必要な情報を行単位で出力します。

summayの結果を見てみる

rdsrSummary.csvファイルの出力結果は次のようになります。

最初に基本的な情報が表示されます。

以下ではCTのX線量に関する情報が単位とともに表示されています。

おわりに

今回、2022年4月24日に日本放射線技術学会北海道支部第78回春季学術大会　医療情報セミナーで発表したものを掲載しました。

DICOM SRの読み込みはPydicomのサンプルrtplan.srで試していましたが、今回RDSRは初めての経験でした。

読み込むまでは良かったのですが、階層が6階層まで下りて行く必要があり、Pydicomの標準機能ではうまく情報を取り込むことができませんでした。

CTのほかにAG（血管造影）のRDSRも解析しましたが、同じ結果でした。

本来はDICOM規格をしっかり理解したうえでプログラミングすべきでしたが、短期間に結果を出す必要があり、力ずくの対応でした。

したがってまだ気づかぬバグなどが潜んでいるかもしれません。プログラムの使用には注意をお願いします。

登録: 投稿 (Atom)

マイブログ リスト

医療言語処理講座

2022年10月19日水曜日

Pyゼミ0.06 備忘録 Ubuntu20.04にCUDA11.4、TensorFlowをインストールする

はじめに

1）Ubuntu20.04インストール後 ソフトウェアアップデート

2)CUDA＆Driverインストール

3)CUDAパス設定

4)cuDNNのインストール

5）Anacondaのインストール

6）TnesorFlow と Kerasのインストール

7)OpenCVのインストール

8)Pydicomインストール

9)DCMTKのインストール

10)ImageJのインストール

11）その他の設定

○メニューに空のドキュメントを追加する

○エディタ

おまけ）Dropboxのインストール

2022年5月6日金曜日

Pyゼミ1.08 PET画像を読む

医療画像の画素値を定量評価に使う

プログラミング スタート

PET画像を読み16ビットの画素値を取り出す

実行してみましょう

結果画像を見てみましょう。

apply_modality_lutを使ってみる

実行してみましょう

8ビットでウインドニングしてしまう

おわりに

2022年5月2日月曜日

Pyゼミ1.0 番外編 DICOM RDSRを読む

Pydicomを使ってRDSRを読んで必要な情報を取り出してみる

DICOM RDSRの内容を見てみる

少し解析のフローをまとめると

RDSRを読むPythonプログラム

データエレメントの深さ求めるメソッドとデータエレメントの各要素求めるメソッドについて

１．データエレメントの深さを求めるgetDepthメソッド

メインのプログラムを読む

実行してみる

summayの結果を見てみる

おわりに

マイブログリスト

Pyゼミ0.06　備忘録　Ubuntu20.04にCUDA11.4、TensorFlowをインストールする　

1）Ubuntu20.04インストール後　ソフトウェアアップデート

Pyゼミ1.08　PET画像を読む

プログラミング　スタート

Pyゼミ1.0　番外編　DICOM RDSRを読む